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Model granular assemblies, in which grains are assumed rigid and frictionless, at equilibrium under 
some prescribed external load, are shown to possess, under generic conditions, several remarkable 
mechanical properties, related to isostaticity and potential energy minimization. Isostaticity - the 
uniqueness of the contact forces, once the list of contacts is known- is established in a quite general 
context, and the important distinction between isostatic problems under given external loads and 
isostatic (rigid) structures is presented. Complete rigidity is only guaranteed, on stability grounds, in 
the case of spherical cohesionless grains. Otherwise, the network of contacts might deform elastically 
in response to small load increments, even though grains are perfectly rigid. In general, one gets an 
upper bound on the contact coordination number. The approximation of small displacements, that 
is introduced and discussed, allows to draw analogies with other model systems studied in statistical 
mechanics, such as minimum paths on a lattice. It also entails the uniqueness of the equilibrium state 
(the list of contacts itself is geometrically determined) for cohesionless grains, and thus the absence 
of plastic dissipation in rearrangements of the network of contacts. Plasticity and hysteresis are 
related to the lack of such uniqueness, which can be traced back, apart from intergranular friction, 
to non-reversible rearrangements of small but finite extent, in which the system jumps between two 
distinct potential energy minima in configuration space, or to bounded tensile forces, deriving from 
a non-convex potential, in the contacts. Properties of response functions to load increments are 
discussed. On the basis of past numerical studies, it is argued that, provided the approximation of 
small displacements is valid, displacements due to the rearrangements of the rigid grains in response 
to small load increments, once averaged on the macroscopic scale, are solutions to elliptic boundary 
value problems (similar to the Stokes problem for viscous incompressible flow). 

PACS numbers: 46.10.+z,05.40.+j,83.70.Fn 



I. INTRODUCTION 
A. Motivations 

A large research effort, both in the statistical physics 
and the mechanics and civil engineering communities, is 
currently being devoted to granular materials, aiming in 
particular at a better understanding of the relationships 
between grain-level micromechanics (intergranular con- 
tact laws) and macroscopic behaviours (global equilib- 
rium conditions, constitutive relations) P, @, H[ 

This aim -the traditional program of Statistical Me- 
chanics - is far from fully achieved in dense granular 
systems near equilibrium, for one is facing at least two 
fundamental difficulties. 

Firstly, the non-smooth character of contact laws, that 
involve unilaterality and, possibly, dry friction, is a com- 
mon feature of granular assemblies that endows them 
with a high level of disorder and a high sensitivity to 
perturbations. Tiny motions might significantly affect 
the way forces are transmitted, since contacts between 
neighbouring grains might open or close (and the slid- 
ing or non-sliding status of closed ones might change). 
Hence the characteristically heterogeneous aspect of force 
transport in dense granulates: large forces are carried by 
a network of preferred paths (the "force chains") while 
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some grains or sets of grains carry but vanishing ef- 
forts ("arching effect"). The histogram of contact forces 
spans a wide range. These phenomena have been ex- 
perimentally observed thanks to techniques like photoc- 
lastic stress visualization 0, H[ and carbon paper print 
analysis [|| 0] • They have also been studied in numeri- 
cal simulations [1, Q , and some attempts of theoretical 
descriptions have been proposed [lo| . Such peculiar as- 
pects of granular systems render more difficult the ref- 
erence to existing models from other fields. Indeed, a 
recent trend in the physics literature on static granular 
systems [IJ EH EH E3 insists on their difference with 
ordinary, elastic solids, and suggests, instead of resort- 
ing to macroscopic displacement or strain variables, to 
search for direct relations betwen the components of the 
stress tensor. 

The second basic difficulty stems from the incomplete 
knowledge of the mechanical properties of granular sys- 
tems, especially those ruling the dynamics. When a gran- 
ular sample is submitted to some prescribed external ac- 
tions that are sufficiently slowly changing in time, its 
evolution is customarily described as an ordered set of 
equilibrium states that are successively reached, with lit- 
tle or no dependence on physical time. The physical pro- 
cesses by which kinetic energy is dissipated are, however, 
most often somewhat mysterious or poorly characterised. 
They are, in the framework of the quasi-static description 
we have just mentioned, implicitly regarded as irrelevant. 
One might wish to assess the validity of such an assump- 
tion. Numerical simulations, that have to adopt some 
rule to move the grains, could in principle allow useful 
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investigations of the influence of the dynamics. However, 
in view of the practical difficulty to obtain representa- 
tive configurations close enough to equilibrium within a 
reasonable computation time, they sometimes resort to 
non-physical parameters, and pick up the dynamical rule 
among the restricted range of those that allow tractable 
calculations. 

This paper addresses both those basic concerns, in the 
following way. Simplifying assumptions are introduced 
(we consider, e.g., rigid frictionless grains), thus restrict- 
ing our attention to a certain class of model systems, 
that are however argued to exhibit the same qualitative 
behaviours as more realistic ones. Those systems are suit- 
able candidates to test, most easily by numerical means, 
some recently proposed models and speculations, at the 
expense of rather extensive numerical computations. The 
purpose of the present article is not, however, to present 
new results of numerical simulations. We shall state and 
establish, rather, with a fair level of generality, some 
basic properties of such systems, and study their qual- 
itative consequences in terms of macroscopic mechanical 
behaviour. This analysis will shed some light on some 
analogies and differences with other previously studied 
problems in statistical mechanics, such as directed 'poly- 
mers' in random environments and percolation models. 
It will also, along with the exploitation of past numerical 
results on a simplified model [1, 0j| Ell US EH > allow us to 
investigate the possible origins of some macroscopic fea- 
tures of granular mechanics, that are classically modelled 
with elastoplastic constitutive laws 13. [191. a nd to discuss 
other recently proposed approaches [ill Il2l [13 , G3 . 

We will show that mechanics is to a large extent de- 
termined by geometrical aspects (steric exclusion), thus 
partially answering concerns about the role of dynamical 
parameters. Finally we will discuss the status of dis- 
placement and strain variables in quasi-static assemblies 
of rigid grains, and give perpectives for future investiga- 
tions. 



B. Synopsis. 

The paper is composed of two main parts. 

First, sections II to V introduce useful definitions and 
state basic properties that are necessary for the deriva- 
tion of the main results. 

Thus, section II presents useful definitions and me- 
chanical properties of static granular systems, i.e., col- 
lections of rigid bodies essentially interacting via point 
forces mutually exerted on their surfaces. Those notions, 
that include the theorem of virtual power, generalized 
forces and velocities for collective degrees of freedom, 
and the degree of indeterminacy of forces and of veloc- 
ities, are not always familiar in the condensed matter 
physics community. Section III introduces the potential 
energy minimization problems for various simple friction- 
less contact laws. Section IV defines the approximation 
of small displacements, a modelling step of both technical 



and conceptual importance, since it allows, in particular, 
an analogy with problems of scalar transport on discrete 
networks, as explained in section V. 

Once those essential ingredients made available, the 
second part of the paper (sections VI to IX) establishes 
the main results and discusses their consequences, with 
reference to previous theoretical and numerical work, and 
to known aspects of the mechanical behaviour of granular 
materials. 

Section VI is devoted to the generic isostaticity prop- 
erty of equilibrium states in systems of rigid grains that 
may only exert normal contact forces on one another. We 
then prove and discuss (section VII) the uniqueness of 
the equilibrium state in cohesionless systems within the 
approximation of small displacements, and compare the 
determinatin of equilibrium states of such systems with 
other mechanical or scalar transport problems. Section 
VIII introduces the additional requirement of stability, 
outside the approximation, which is dealt with, in the 
absence of friction, in terms of potential energy minimiza- 
tion. In some restricted models, this allows to conclude 
to the isostaticity of the structure, a stronger property 
than mere isostaicity of the problem under a given load. 
It is then possible to discuss the possible origins of plas- 
ticity in systems of frictionless grains and the form of the 
mechanical response to small load increments. The paper 
ends with concluding remarks (X) on the role of displace- 
ments and strains in granular materials and suggestions 
for future research. 



II. BASIC DEFINITIONS AND PROPERTIES. 



We are interested in the modelling of large packings 
of solid bodies (grains), in equilibrium under some pre- 
scribed external forces. Grains are assumed to interact 
via point forces mutually exerted on their surfaces, which 
means that the distribution of stress on their areas of 
contact or of influence can effectively be viewed as local- 
ized at a point, on the scale of the whole grain. Apart 
from this reservation, that excludes flat or conforming 
surfaces [56| . grains might have arbitrary shapes, and 
our considerations apply to spatial dimension d equal to 
2 or 3, although most examples will be taken with two- 
dimensional systems of discs. Note that we do not require 
interacting grains to touch one another at this stage. We 
mostly restrict our attention here to frictionless bodies, 
i.e., such that contact forces are normal to the grain sur- 
faces. This might look like a severe limitation, but wc 
shall argue that such simplified systems do possess the 
generic properties of granular media. We shall also as- 
sume, unless otherwise specified, that the grains behave 
as rigid undeformable objects. 
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A. System, external forces 

Wc consider a set of n grains, labelled with indices i, 
with 1 < i < n. In each of them we arbitrarily choose 
a 'center', which might e.g., coincide with its center of 
mass. In the case of spherical grains it is of course con- 
venient to take the geometrical center of the sphere. The 
(d-dimensional) velocities of those centers, (Vj)i<,-< n , to- 
gether with the d'-dimensional (with d! — £^£L_L!) angular 
velocities {£li)\<i<m make up the kinematic degrees of 
freedom of the whole system, thus labelled by couples of 
indices (i, a), with 1 < i < n and 1 < a < d+d' = , 
We denote as / the set of such couples. If a > d, Vi iCt is 
now a notation for f2j iQ _,2. Boundary conditions are of- 
ten enforced by prescribing the motion, or the absence of 
motion, of walls. Those might be regarded as solid bod- 
ies, or particular 'grains' themselves. In the following we 
shall sometimes write down large 'velocity vectors' that 
gather all Nf kinematic degrees of freedom of the system, 
then denoted, with a single index, as (w M )i< M <Ar / . 

It might also be convenient to keep some grain coordi- 
nates fixed (thus choosing one particular Galilean frame) , 
i.e., to impose, for all couples i, a belonging to some sub- 
set Iq of /, Vi. a = 0. Indices /i are then renumbered, and 
Nf is reduced accordingly, to label and to count the free 
kinematic parameters. Another classical way to impose 
some boundary conditions is to require, for all i, a in 
some subset I\ of /, Vi >a to depend linearly on one or 
several parameters, e.g.: 

(i,a)€li Vi, a = A iiCC Xi, (1) 

introducing some collective 'generalized velocity' Ai. 
Once again, in such a case, Nf is reduced to count el- 
ements of I\ (Iq U Ji), plus Ai. 

At least locally, it is possible to regard velocities and 
generalized kinematic parameters (like Ai in eqn. [I} as 
time derivatives of spatial coordinates, which we shall do 

in the following, thus writing, e.g., Vi a = — As we 

at 

are only interested in those properties that do not de- 
pend on dynamics, grain trajectories might as well be 
described by any parameter, not necessarily by physi- 
cal time. In the case of kinematic constraints of type 
[U parameters A ijQ will be regarded as fixed, although 
positions of the grains and the walls change. One then 
defines a generalized coordinate Ai, such that = Ai. 
Just like for velocities, the compact notation (x M )i< M <Ar / 
refers to the whole set of positional coordinates. 

External forces and torques may at will be exerted on 
the grains that are free of kinematic constraints. We shall 
use the same notations as for velocities, writing down 
large iV/-vectors of 'external forces' (some of their coor- 
dinates standing, actually, for torques), as {Ffi Xt )i<^<N r 
At equilibrium, they are of course to be balanced by in- 
ternal forces (^"'H^jv,: 

(1 < n < N f ) F* xt + Fl nt = 0. (2) 



In order to enforce constraints of type [IJ some external 
efforts have to be exerted on the concerned bodies. On 
requiring the power of such efforts to be balanced by 
that of internal forces {F™ t )i<ij,<.Nf, onc identifies the 
generalized force conjugate to Ai as 

(i,a)eli 

We just used the power to find generalized forces: this is a 
manifestation of the duality between forces and displace- 
ments or velocities, which will be repeatedly exploited in 
the sequel. The iV/-dimensional vector space T of ex- 
ternal forces, is, by construction, to be regarded as the 
dual space, in the ordinary sense of linear algebra, of the 
Ay-dimensional space V of kinematic degrees of freedom. 

In general, it should be appreciated that the appro- 
priate mathematical description of configuration space is 
not ]R Nf with its Euclidean structure, but, due to rota- 
tional degrees of freedom, an iV/-dimensional manifold, 
on which (a; A ,)i< M <Ar / is a set of (local) curvilinear coor- 
dinates. V and T are respectively the tangent and cotan- 
gent vector space at a given point, and depend on that 
point. Thus the definition of 'constant velocities', or of 
'constant forces' requires some care. However, these dif- 
ficulties are inessential in our subsequent treatment, and 
we shall assume 'constant external forces' are applied, 
and derive from a potential energy: 

W = -J2 F ™%- ( 4 ) 

It is easily checked that such a definition is devoid of 
ambiguity in the following important cases. 

• The set of grain center positions, as opposed to 
grain orientations, define a 'flat' space, on which 
constant vectors and covectors are unambiguous. 
Whenever external efforts are not sensible to ori- 
entational coordinates, as in the case of gravity (if 
the grain 'centers' are their centers of mass), one 
may therefore 'apply constant forces'. 

• Anticipating on part IV, the approximation of small 
displacements assumes that the manifold might lo- 
cally be replaced by its flat tangent space. 

The complete JV/-vector of external forces is referred 
to as the load. Sometimes, it is convenient to deal with 
parametrized sets of loads. When the direction of the 
load is fixed, while its intensity might vary, one has a one- 
parameter loading mode. In such a situation, all external 
force components are kept proportional to a single load- 
ing parameter Q, and a generalized velocity conjugate to 
Q, A, can be identified on equating the power of external 
forces with the product QX. A is some linear combination 
of the kinematic degrees of freedom (t | M )i< A ,<Ar / , and the 
time derivative of a generalized coordinate A, equal to 



4 



the same combination of coordinates {%ij,)i<(i<Nf The 
potential energy is then simply 



W = -QA. 



(5) 



Let us now illustrate those notions with simple exam- 
ples, that will be repeatedly used in the following. Sys- 
tems A and B are packings of discs that are placed on 
the sites of a regular triangular lattice. (Later on, we 
shall allow for a slight polydispersity of the grains. They 
might move, gain or lose contacts with their neighbours, 
and the lattice might be slightly distorted). System A 
(figtU) is a pile with slope inclined at 60 degrees with re- 
spect to the horizontal direction. Each disc is submitted 




FIG. 1: System A : a pile under gravity. The bottom bound- 
ary conditions are explained in the text 



to its own weight, except those of the bottom row, which 
collectively set the boundary condition. One might keep 
them fixed at regularly spaced positions, imposing, say 
(numbering them as on the figure, and denoting as a the 
lattice spacing) 



(1< i < 



Xi =(i - 1)(1 - Ai)a 
Vi =0 



(6) 



allowing for a horizontal deformation parameter Ai. One 
may also require them to stay on the horizontal axis y = 
and satisfy 



-Ai(i - l)a, 



(7) 



with a free kinematic parameter Ai. According to eqn.[3J 
the generalized force conjugate to Ai is 



(8) 



These two slightly different boundary conditions (BC) 
are respectively abbreviated as BC1 and BC2 in the fol- 
lowing. 

System B (fig. [5J is a hexagonal sample of the same 
material. It is submitted to external forces on the pe- 
riphery, which mimic hydrostatic pressure. 




FIG. 2: System B : a hexagonal sample. Arrows depict ex- 
ternal forces applied on peripheral discs. 

System C (fig-EJ) is a disordered collection of discs with 
a larger polydispersity. It is embedded within a circular 
wall the radius R of which might change. One controls 
the generalized force conjugate to Ai = viz. 



Qi — E & 



(9) 



where the sum runs over all particles i exerting forces 
normally onto the wall. 



B. The structure: a set of bonds. 

The definitions we introduce here pertain to one spe- 
cific configuration of the grains, with the positions and 
orientations fixed. 

We call 'bonds' the pairs of neighbouring grains that 
may exert a force on one another. We require this force 
to be concentrated at the point of each grain which is 
the closest to the other one, and directed normally to 
the surface. [5 7I| The more general case of arbitrary bond 
forces will be briefly evoked later. Note that we neither 
require the grains that are joined by a bond to be in 
contact, nor impose any sign constraint on the force. We 
thus define, somewhat arbitrarily at this stage, N such 
bonds as depicted on figure [4] alternatively labelled with 
an index I, 1 < I < N, or with the pair of labels of the 
two grains they join. If bond I connects i and j, n/ or riy 
denotes the unit vector that points from i to j, normally 
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FIG. 3: System C : a disordered packing surrounded by a cir- 
cular wall that might uniformly expand or shrink, as indicated 
by the small arrows. 




FIG. 4: Two grains i and j joined by a bond, hij is the 
minimum distance between their surfaces, measured where a 
common normal unit vector is ny. Vector Ry (respectively 
Rji) points from the center of i (resp. of j) to the point of 
its surface that is closest to j (resp. to i) 

to the surfaces of both grains where the distance between 
them , hij, is the smallest. Ry is the vector joining the 
center of grain i (origin) , to the point on its surface that 
is closest to grain j (extremity) . This contact zone might 
transmit a normal force, along n.y , of magnitude /y that 
will be counted positively when the grains repell each 
other. Once this set of bonds is defined, it is referred to as 
the structure. The set of bonds defined by intergranular 
contacts (hij — 0) will be called the contact structure. 
As a consequence of the definition of a structure, 



the form of internal forces ((F- nt )i<i<„) and torques 
((r™*)i<i<„) in the system is specified: they linearly 
depend on bond forces /y, as 

1 . 2^ J.j " .j A Hij 

Given the load (F^ xt )i<^<N f , equilibrium requires, in 
view of eqns. [TUl and 121 that the bond forces (/;)i<;<jv 
satisfy equations of the form 

N 

defining a linear operator, H : JR — > T. Bond forces 
{fi)i<i<N then said to be statically admissible with 
the load {F^ xt )i<ti<N f - Bond forces that are statically 
admissible with a load equal to zero (in equilibrium with- 
out any external action) are the elements of a subspace 
So of PR^, the null space of operator H . Its dimension, 
that we denote as h, is the number of linearly indepen- 
dent such self-balanced sets of internal forces, or, in other 
words, the degree of indeterminacy of bond forces in the 
system (also called the degree of hyperstaticity) . If not 
empty, the set of statically admissible bond forces is an 
affine space of dimension h. 

The relative normal velocity of the grains i and j joined 
by a bond is 

SVij = n, , • (V, - V, + fli A Rij - Slj A R J2 ) , (12) 

with the convention that it is positive when the parti- 
cles are approaching each other. Eqn. 1121 defines a linear 
operator, G, acting on V into IR . The range of G is 
the subspace C of compatible relative normal velocities, 
i.e., those N- vectors for which one can effectively find 
values for the velocities, relations [T2l being satisfied. The 
null space of G is the vector space M of 'mechanisms', 
also called 'floppy modes', i.e., motions that do not alter 
the lengths hi of the bonds. Its dimension, denoted as 
k in the sequel, is the number of independent such mo- 
tions, or, in other words, regarding the bonds as rigid, 
the degree of indeterminacy of velocities, also called de- 
gree of hypostaticity. Imposing the condition SVij = 
in all bonds of the structure restricts the possible values 
of velocities (v fJi )i< fl <N f to a vector space of dimension 
k. Depending on the type of load and boundary con- 
ditions, the whole set of grains might keep some overall 
rigid body kinematic degrees of freedom. System B, for 
instance, has 3 independent such motions, as any solid 
body in 2D. If fco < d(d + l)/2 denotes the number of 
such particular motions allowed by the boundary condi- 
tions, the system is said to be rigid when it does not have 
other mechanisms, i.e., when k = fco. 

An important and useful result, the classical theorem 
of virtual power states the following. Let (SVi)kkn be 
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any element of C, corresponding to the velocity vector 
( v n)i<ti<Nf, and let (fi)i<i<N be a set of bond forces 
statically admissible with the load (F^ xt )i<n<N- One 
then has: 

N N f 

J2fiSVi = J2K xt ^- ( 13 ) 

Equalitv ll31 for arbitrary ('virtual') equilibrium set of in- 
ternal forces and velocities, stresses the geometric mean- 
ing of forces and the mechanical meaning of velocities. 
It is easily established in two steps: first use the force 
balance equations in the right-hand side; then transform 
the sum over degrees of freedom into a sum over bonds. 

As a direct consequence of the theorem, one deduces 
that operator H is in fact (as one might check directly, 
reading the matrix elements in eqns.[T2"land llip the trans- 
pose of G : H = G T . This follows from the sequence of 
equalities 

(f\SV) = (f\Gv) = (Hf\v) = (G T f\v) , 

valid for arbitrary v (such that Gv = Sv) and / (such 
that Hf — F ext ), in which a bracket notation is used for 
scalar products. Consequently, So, the null space of G T , 
is the orthogonal complementary to C, the range of G, in 

So = C x . (14) 

Thus to check that some values SVi that one might try to 
assign to the relative normal velocities are compatible, 
it is sufficient to ensure the orthogonality of iV-vector 
(8Vi)i<i<n to all TV-vectors of self-balanced bond forces 
(or a spanning subset thereof): 

TOi<z<jv -L S . (15) 

One thus uses forces (elements of So) as cofactors in a 
set of geometric compatibility conditions. 

Recalling k (the number of mechanisms) is the dimen- 
sion of the null space M of G, one has 

Nf =fc + dim(C). 

As h = dim(So), from 1141 one also has: 

N = /i + dim(C). 

Elimination of the dimension of C from those two equal- 
ities yields the following relationship between the degree 
of hypostaticity, k, the degree of hyperstaticity, h, the 
number of bonds, N, and the number of degrees of free- 
dom, Nf. 

N + k = N f + h. (16) 

As we will check on examples below, relation [TH] holds 
whatever the choice of the list of bonds between objects, 
although it is of course desirable in practice to define 



bonds according to the interaction law. One may, for 
example, declare a bond to join two grains whenever their 
surfaces are separated by a minimum distance smaller 
than some threshold ho > 0. The choice of a larger ho, 
thereby increasing AT, will decrease k and/or increase the 
degree of hyperstaticity h. 

Let us remark that the properties we have just dealt 
with in the case of bonds that carry normal forces, are 
very easily generalized to the case of arbitrary contact 
forces, at the cost of minor modifications. Relative nor- 
mal velocities and normal contact forces are replaced by 
d- vectors, IR^ replaces IR^, equalities IT51 (with, now, 
a scalar product within the sum in the left-hand side) 
and [TJ] are still satisfied. Instead offTBl one ends up with 
dN + k = Nf + h. Adding friction increases h and/or 
decreases k. 

Returning to frictionless systems, the case of spheres or 
discs deserves a special treatment: no normal force is able 
to exert any torque, and all rotational degrees of freedom 
are therefore mechanisms. It is convenient to ignore them 
altogether. Their number n d ( rf ~ 1 ) ( n j s t ne number of 
particles) is then subtracted both from Nf and from k, 
and eqn llGI still holds. Such granular systems are then 
analogous to 'central-force networks': networks of freely 
articulated bars, or systems of threads tied together, in 
which only the translational degrees of freedom of the 
nodes matter. One should be aware, however, that the 
presence of friction reinstates rotations into the problem. 

We now illustrate the notions and properties intro- 
duced in this section with examples of structures defined 
in systems A, B and C, ignoring, as explained just above, 
disc rotations. 

First consider system B. Three different structures are 
apparent on figure O The first one, that we denote as 
SB1, is the set of bonds that are drawn as thick lines; the 
second, SB2, contains all bonds of SB1, plus those that 
are drawn with thin continuous lines on the figure; and, 
finally, the third structure, SB3, comprises all possible 
bonds between nearest neighbours in the system, i.e., all 
those of SB2 plus the dotted lines. Ignoring rotations, 
one has Nf = 2n = 38. 

Structure SB3 is a set of rigid triangles sharing com- 
mon edges with their neighbours. It is devoid of mecha- 
nisms, except the 3 overall rigid body degrees of freedom 
of the system. Thus k = 3. N = 42 bonds are present. 
In view of eqn. 1161 one has h — 7. One can exhibit 7 lin- 
early independent systems of self-balanced normal forces, 
as follows. The small structure, with 12 bonds, involving 
7 discs, depicted on fig. allows to define one such set of 
forces. Noting that 7 such patterns are present on SB3 
(centered on discs 5, 6, 9, 10, 11, 14 and 15), the right 
count is reached. 

Structure SB2 is made of N — 35 bonds. It can be 
shown (on studying the properties of the corresponding 
matrix G) to be devoid of self-balanced sets of forces, 
h = 0, and of mechanisms other than rigid body motions, 
k = 3. Thus N + k = N f + h. 

Structure SB1, comprising N — 25 bonds only, still 
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FIG. 5: A set of self-balanced normal forces. The 6 bonds 
of the regular hexagon perimeter carry some normal force /, 
while the 6 ones involving the central disc labelled 1 carry the 
opposite force. 

has h — 0. According to eqn. [IHJ it should possess 10 
additional independent mechanisms. 2 of them are due 
to disc 10, which is now completely free. 4 others involve 
discs 5, 9, 12, and 19, which are still free to move in one 
direction. In the case of a divalent disc like 5, this is due 
to the exact alignment, on the regular lattice, of bonds 
4-5 and 5-6. Four less trivial mechanisms are more collec- 
tive. One of them is shown on figure [5J Two structures, 




FIG. 6: A collective mechanism on structure SB1. Arrows 
represent disc velocities. 

SA1 and SA2, are denned, on fig.[TJ in system A. SA1 is 
made of all bonds drawn with continuous lines, and SA2 



contains, in addition, the two bonds drawn with dotted 
lines (19-24 and 32-34). Depending on the boundary con- 
dition, discs 1 to 8 either possess collectively one degree 
of freedom (for BC2) and then Nf — 57, or none (for 
BC1) and N f = 56. 

SA2 has 57 bonds. It is devoid of mechanism (k = 0) 
for whatever BC. For BC2, one also has h = and eqnfTrjl 
holds as an equality between the number of bonds and 
the number of degrees of freedom. For BC1, one has 
h = 1. Indeed, one may recognize, in the bottom left 
corner of the pile, with discs 1, 2, 3, 9 and 10, part of 
the hyperstatic pattern of fig. O With BC1, one needs 
not care about equilibrium of discs 1, 2 and 3 that are 
perfectly fixed. A system of self-balanced bond forces is 
thus found on attributing a common value to the normal 
forces in bonds 1-9, 9-10, 10-3, and the opposite value 
to the normal forces in bonds 2-9 and 2-10. In the case 
of BC2, those forces do not balance, since the equilib- 
rium equation for the collective degree of freedom of the 
bottom row (a combination of eqns. [8landfT0|) is not sat- 
isfied. As to SA1, it has the same properties as SA2, 
with 2 additional mechanisms (collective ones like that 
of fig. ED. 

Consider now structure SC that is shown, in system 
C, on figure [3J with the lines connecting disc centers, 
or joining discs to the wall, that define N = 70 bonds. 
Taking into account the degree of freedom of the wall, one 
has Nf = 2n + 1 = 75. One may show h = 0. Thus one 
has k = 5. Two discs (10 and 14) are entirely free, hence 
4 mechanisms. The missing one is a global rotation, as a 
solid body, of the set of all particles around the center of 
the circular container, the wall remaining immobile. Such 
a motion would not be possible if the same boundary 
condition was used with another container shape. 



C. The problem: the structure and the load. 

Once a list of bonds is chosen, thus defining the struc- 
ture, we shall refer to the situation of the structure sub- 
mitted to a given load as 'the problem'. 

Solving the problem would mean finding the motion or 
equilibrium state of the system (determining, e.g., new 
equilibrium positions and intergranular forces), once the 
load, from an initial state of rest with no external force, 
has been applied. We are not, of course, able to do that 
at this stage, since no contact law relating the forces to 
the relative motion of neighbouring particles has been 
introduced. The only information available is that the 
internal forces are required to belong to some vector space 
that is known once the structure is defined, and to be 
exerted on given points on the grain surfaces. 

It is said that the load is supported by the structure 
if its application leads to an equilibrium state in which 
internal forces, carried by the bonds of the structure, 
balance the external ones. 

We can state a necessary condition for the load to be 
supported: it must be possible to find statically admis- 
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sible intergranular forces. Necessarily, the 7V/-vector of 
external forces must lie in the range of operator G T , i.e., 
it must be orthogonal to the null space M of G: 

(F^iiK^N,) ± M. (17) 

This simply means that if the load is to be supported, it 
must not set the mechanisms into motion. Such a load 
is said to be supportable. All supportable loads are not 
always supported. 

By definition, the backbone of a structure is the set of 
bonds Iq such that a list of statically admissible internal 
forces (fi)i<i<N exists with fi ^ 0. In the following we 
shall also refer, as 'the backbone', to the set of grains 
reached by such bonds. 

In general, a full mechanical characterization of the 
equilibrium properties of the system requires some con- 
stitutive law in the contacts. However, there are inter- 
esting situations in which 

• condition [T7l being fulfilled, the load is supportable; 

• if it is supported, then all intergranular forces are 
uniquely determined by the equations of equilib- 
rium. 

These two conditions define an isostatic problem. 

Further restrictions on internal forces are often en- 
forced in the form of inequalities. The definition of a 
supportable load is then modified accordingly, impos- 
ing additional conditions, to be satisfied simultaneously 
with 1 171 Their consequences will be discussed in sections 
III and VII. 



D. Isostaticity: various definitions. 

In section VI we shall see that equilibrium configura- 
tions of assemblies of rigid frictionless grains interacting 
via contact forces only are generally such that the prob- 
lem is isostatic. 

Here, we first insist on the difference between an iso- 
static problem, as defined just above, and an isostatic 
structure, to be defined below. Once condition [T7] is sat- 
isfied, the set of possible bond forces is an affine space 
of dimension h. One has an isostatic problem if both 
conditions 1171 and h — are fulfilled. Some mechanisms 
might still exist in the structure (k ^ 0), provided they 
are orthogonal to the load direction. 

Structure SAl (figure [1]), with discs exactly centered 
on the sites of a regular triangular lattice, is such that 
the problem, denoted as PA1 in the following, defined 
with BC2 and the following load: [58| 



(9 < i < 36) Ff 



-pe y 



(18) 



where p is the weight of one disc and e y is the vertical 
upwards unit vector, is isostatic, although 2 mechanisms 
are present. 



Analogously, structure SB1, along with the load shown 
on figure [21 defines an isostatic problem PB1 in spite of 
the k = 10 mechanisms. In particular, the load direction 
(provided discs sit right on the regular lattice sites) is 
exactly orthogonal to the velocity vector represented on 
figure [H Structure SC, submitted to the following load: 



(1<*<37) 



= 

0?, 



(19) 



where a prescribed value Q\ is imposed to generalized 
force Q\ defined in eqn. [9l yields an isostatic problem. 

Isostatic structures, on the other hand, are such that 
all problems are isostatic, whatever the choice of the load. 
More precisely, one requires all loads orthogonal to the 
overall rigid-body degrees of freedom to be supportable 
with a unique determination of internal forces. Equiva- 
lently, both conditions h = and k = fco are to be satis- 
fied. Both the degree of hyperstaticity and the degree of 
hypostaticity (excluding rigid-body motions) should be 
equal to zero. This entails the well-known condition 



N = Nf — fc , 



(20) 



stating that the number of equilibrium equations (Nf — 
fco) is equal to the number of unknowns (N). 

Equality [201 is a necessary condition for the structure 
to be isostatic, not a sufficient one. For example, in 
the structure defined by the addition of the bond joining 
discs 19 and 24 to SAl with the first boundary condition 
(BC1), one has fc = 0, N = N f = 56, while h = k = l. 

Structure SA2, with BC2, is isostatic. SB2, with 
Nf = 38 and fco = 3, is isostatic. As to SC, it would 
be isostatic upon removal of grains 10 and 14, only if the 
global rotation of the set of grains with respect to the 
wall were ignored. Of course, all those structures, as we 
are dealing with discs, are only isostatic if rotations are 
ignored. Only problems with no external torque exerted 
on the grains are isostatic. This should be remembered 
on comparing h and fc with and without friction in such 
systems. 

As we shall see, isostatic problems, rather than iso- 
static structures, naturally occur in some model granu- 
lar systems. The distinction is relevant, for it accounts 
for disconnected or 'dangling' parts in disordered struc- 
tures like SC, and for the peculiarities of lattice models. 
Moreover, some systems can also spontaneously, as we 
shall see, select a non-rigid (fc > fco) equilibrium config- 
uration. 



E. Generic versus geometric properties. 

The distinction between isostatic problems and iso- 
static structures should not be confused with another 
one: that between geometric and generic isostaticity. We 
have used a geometric definition of a structure, as asso- 
ciated to one particular position of the system in con- 
figuration space, and accordingly the definition we gave 
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is that of geometric isostaticity. A topological one can 
be introduced which, irrespective of particle positions, 
is only sensitive to the connectivity of the network of 
bonds. In the case of spheres or discs, when rotations 
can be ignored, this amounts to regarding the structure 
as a graph: a set of edges (bonds) joining at vertices 
(grains). Operator G, spaces Sq, M and their dimen- 
sions h and k smoothly depend on the coordinates of the 
grains, via vectors riy and Ry . However, the rank of a 
parameter-dependent matrix stays at its maximum ex- 
cept for special values of the parameters. Equivalcntly, 
the dimension of the null space is generically equal to its 
minimum value. Applying this to both G and G T , one 
may define the generic degree of indeterminacy of veloc- 
ities (with due account to the kg rigid-body degrees of 
freedom) k and the generic degree of indeterminacy of 
forces h as the respective generic (minimum) dimensions 
of their null spaces. This allows to define a suitable iso- 
staticity notion for topological structures: a generically 
isostatic structure is one for which both numbers h and 
k — ko are equal to zero. 

It follows from the definitions that a geometrically iso- 
static structure is always, once regarded as a topolog- 
ical structure, a generically isostatic one, but that the 
reciprocal property is not true. Ref. [2(| gives a coun- 
terexample for a system of discs (like systems A and B, 
equivalent to a network of articulated bars) on the regu- 
lar triangular lattice. In specific configurations (like that 
of a regular triangular lattice), one might exceptionnally 
have h = k — fco > on generically isostatic structures. 

In two dimensions, there exists some powerful algo- 
rithms [U to evaluate the generic degrees of force 
and velocity indeterminacy in central-force networks (or 
systems of frictionless discs). Such computational meth- 
ods only deal with connectivity properties, they do not 
manipulate floating-point numbers and are therefore de- 
void of numerical round-off errors. They were success- 
fully applied to systems of up to 10 6 nodes. However 
they are of course unable to compute position-dependent 
quantities like force values. 



III. CONTACT LAW AND POTENTIAL 
MINIMIZATION. 

So far, the only restriction on intergranular forces was 
that they should be normal to the grain surfaces. [59J In 
this section we consider some more specific cases of fric- 
tionless grains, in which some "contact law" , relating nor- 
mal forces to relative positions, is known. This provides 
some limited additional information, that is not sufficient 
in general to predict the grain trajectories once they are 
submitted to external forces, for all dynamical aspects 
are still unknown and the characterization of equilibrium 
might even be incomplete. Our aim is to deduce as much 
as possible on the global properties of the granular assem- 
bly from as little information as possible on the detailed 
mechanical laws of the contacts, in order to stress the 



importance of geometrical aspects. Thus we first present 
the simplest case of rigid, frictionless and cohesionless 
grains, in which contacts simply behave as struts. Then 
we introduce and briefly discuss other possible laws in 
which unilaterality or rigidity constraints are modified or 
relaxed. Most of those frictionless systems possess a po- 
tential energy that is stationary at equilibrium states and 
reaches then a minimum if they are stable. Throughout 
this section, it is assumed that a one-parameter loading 
mode has been defined for varying particle positions and 
orientations, with constant external forces, and that the 
potential energy of external forces, W, can be written in 
the forms of eqns. [4] and [5] 



A. Rigid frictionless grains, no cohesion. 

In this case, the contact law takes the form of the so- 
called Signorini condition: 



if ha > 



fij > 



if hij = 



(21) 



It should be noted that this law does not express a func- 
tional dependence of fi on hi . Let us study the variations 
of W near equilibrium states. First, consider such a state, 
in which some non-negative contact forces /*, in closed 
contacts (hi — 0) balance the external load Q. Let us ap- 
ply the theorem of virtual power with statically admissi- 
ble force set (f*)i<i<N, and arbitrary particle velocities, 
corresponding to relative normal velocities SVi = — ^r- 
and a value A = 4r for the kinematic parameter conju- 
gate to Q. For any I such that /* > 0, the Signorini con- 
dition requires that hi = and one must have 5Vi < to 
comply with the impenetrability constraints. Then, from 



dW 
~~dt 



dA 
dt 



it follows that any motion that does not lead to grain 
interpenetration can only, to first order in t (any param- 
eter on the trajectory in configuration space) increase the 
potential energy. This non-negative first-order variation 
might be equal to zero if SVi = for any active contact I, 
i.e., if a mechanism exists on the backbone of the contact 
structure. Whether the equilibrium state corresponds to 
a minimum of W depends then on the sign of second 
or higher order variations. If the backbone of the con- 
tact structure is rigid, then W is necessarily minimized 
at equilibrium. 

Conversely, let us assume that a configuration of the 
grains has been reached, that locally minimizes W un- 
der the constraints hi > 0. There must then exist some 
non-negative Lagrange multipliers such that, for any 
coordinate x a , 



d 

dx„ 



(QA) 



dhi 
dx a 



(22) 
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Only for such indices I that hi — do the fi take non- 
vanishing values. The partial derivative in the right-hand 
side of [22] is the opposite of matrix element Gz, a , while, 
fromQJ that of the left-hand side is the external force con- 
jugate to x a . Thus, we have just written that parameters 
fi are in fact equilibrium contact forces satisfying[2T] and 
reaction forces stem from geometrical constraints. 

We now introduce a few other related contact laws and 
mechanical models. 



B. Systems with tensile or bilateral forces. 

Networks of rigid strings or cables are analogous to 
frictionless spheres (ignoring their rotations) if the sign of 
forces is reversed and if the distance constraint hi > is 
replaced by hi < 0. The Signorini condition |2"T1 becomes 



fij = if hi 
fij < if hi 



< 

= o, 



(23) 



and the whole treament of the preceding subsection 
straightforwardly applies. 

In the case of non-spherical grains, an analogous sys- 
tem supporting tensile forces is an idealized chain, in 
which 'grain' -chain links- perimeters are free to cross. 
Pairs of neighbouring links (interpenetrating 'grains') ex- 
ert a force on one another, opposing their separation, 
when their intersection reduces to a contact point. 

A bilateral contact law: 



h 



o 



if hij ^ 



fij unknown if h^ = 0, 



(24) 



might model rigid cohesive grains, that 'stick' to one an- 
other. The sticking force might be limited by an unequal- 
ity: 



fij =0 if hij f 
fij > -fa if Kj = 0, 



(25) 



When one simply uses the form [24] assuming the pairs 
that are stuck in contact will not come apart, the conclu- 
sions of subsection III. A still hold, if unilateral conditions 
on relative velocities and displacements are replaced by 
bilateral ones, and if all sign constraints on contact forces 
are removed. Equilibrium configurations are character- 
ized by stationarity of potential energy W. Minimization 
of W ensures stability. A sufficient, but not necessary 
condition for minimization of W is the rigidity of the 
backbone of the contact structure. 

Reciprocally, statically admissible normal contact 
forces naturally appear as Lagrange multipliers associ- 
ated with bilateral constraints ft; = at a potential en- 
ergy minimum. 

However, contact law [55] does not lend itself to a po- 
tential energy formulation. 

Tensegrities [231 ] (with rigid elements) are by definition 
mixed networks of struts (satisfying condition al)) or bars 



(bilateral) on the one hand, and cables (satisfying | 
on the other hand. Their potential energy has the same 
properties as stated above. 



C. Systems with a smooth interaction potential. 

The model of perfectly rigid grains is physically rea- 
sonable when contact deformations (hi < 0) are negligi- 
ble in comparison with any other relevant length in the 
problem. When this is no longer the case, or when one 
wishes to model sound propagation, it is appropriate to 
deal with contact laws that involve elastic deformations, 
e.g., 



■fij 
■fij 





K. 



ij i hij | 



if hij > 
if hij < 0, 



(26) 



in which Kij is a stiffness constant that depends on mate- 
rial properties and on the geometry of contact The 
exponent is m = 3/2 (Hertz law) for smooth surfaces 
in 3D, and other values might model roughness and the 
presence of conical asperities [U HI[ . 

Such contact forces derive from an elastic potential en- 
ergy: 



N 



K, 



w(hi), with w(hi) = 
z — ' m + 1 

1=1 



\h,\ m+1 . 



(27) 



Likewise, rigid cables as introduced in subsection III.B 
could be replaced by elastic ones. That stable equilibrium 
states correspond to minima, in the absence of frictions, 
of the total potential energy 



W + W, 



(28) 



sum of the elastic potentiall27land the potential energy of 
external forces 2] or [5] is an extremely familiar property. 
The Signorini condition might physically be regarded as 
the limit of the interaction law expressed by equation 1261 
when the stiffness constants become very large, or, equiv- 
alently, when the level of intcrgranular forces approaches 
zero. Alternatively, it is mathematically possible to in- 
troduce a regularized contact law of the form [26] as an 
approximation, when contacts are stiff enough, of the 
ideal impenetrability constraint. Such a point of view 
is adopted in optimization theory: the procedure known 
as penalization of the constraints amounts, instead of 
minimizing W subject to impenetrability constraints, to 
searching for unconstrained minima of W + W el . 

Tensile contact forces of limited intensity, as in contact 
law I25[ might result from some attractive interaction of 
finite, but small, range, as depicted on fig. [7] It is in- 
teresting to note that the addition of an attractive tail 
has turned the potential w(h) into a non-convex function 
of interstitial thickness h. At the inflexion point, A, the 
attractive force reaches its maximum f If one pulls, 
with a growing force, on two grains in contact in order to 
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= 0.5 



-0.5 




0.06 



FIG. 7: Interaction potential, with an attractive tail, as a 
function of interstitial thickness h. The curve has an inflex- 
ion point A, corresponding to the maximum attraction force 
(equal to the slope of the dotted line). 



separate them, an instability, in which the contact sud- 
denly breaks open, is reached as the pulling force reaches 
the value /a- When the corresponding intcrgranular dis- 
tance, IfiA-, is so small that it is negligible in comparison 
to all other relevant lengths in the problem, one might 
then replace the smooth attractive potential by contact 
law[25l with /o = /a- On doing so, one loses however 
the possibility to exploit minimization properties. 

We shall see that the potential minimization proper- 
ties have important consequences in terms of the pos- 
sible uniqueness of the equilibrium state under a pre- 
scribed load, and, eventually, as to the possible origins 
of macroscopic plastic dissipation. But, first, we have 
to extend the properties we have stated for velocities 
(or infinitesimal displacements) to small displacements, 
around a given reference configuration. 



able if one wishes to deal with linear problems: adding 
up two displacement fields, for instance, in continuum 
mechanics, is otherwise a meaningless operation. In the 
case of granular systems, it will also lead to a lineariza- 
tion of the problems, for the curvature of configuration 
spaces will be ignored. Its range of validity has to be as- 
sessed a-posteriori, but is of course presumably larger in 
dense systems, where contacts might open and close with 
only tiny changes of the relative positions of neighbouring 
grains. 

Specifically, we assume the coordinates of the grains 
to stay close to reference values. Quantities pertaining 
to the reference configuration will be labelled with a su- 
perscript '0'. It is often convenient, then, to work with a 
fixed structure -the list of contacts that might close, and 
transmit a force, is a-priori known. 

Interstitial thicknesses hi are written as hi — h® — Sui , 
with a relative normal displacement Sui that is linear 
in the grain displacements (and rotations), regarded as 
small quantities. Vectors n^, R^, R^-j are regarded as 
constant, equal to n? , R?-, R^. As they appear as co- 
factors of the displacements, taking their variations into 
account would introduce second order terms. All changes 
of the structure geometry are ignored. Spaces C, 5, oper- 
ators G, G T , are assumed to be the same in the actual as 
in the reference configurations. Displacements are now 
endowed with the same linear algebraic structures as ve- 
locities. G operates on displacements, yielding relative 
normal displacements Sui , the compatibility condition for 
relative normal displacements is the orthogonality to the 
space of self-balanced internal forces So, a theorem of vir- 
tual work can be stated instead of the theorem of virtual 
power, etc. . . 

Within the framework of the ASD, the specificity of 
mechanical problems disappears: as the effect of the dis- 
placements of the grains (variations of the coordinates) 
on the positions (coordinates) themselves are ignored, 
one can find analogies with various other local properties 
of a list of fixed points, nodes or lattice sites. Forces now 
appear as unknown vectors carried by fixed directions, 
and the sum of incoming forces on a node has to vanish. 
Part V introduces the analogy with scalar transport on 
a fixed network. 



IV. THE APPROXIMATION OF SMALL 
DISPLACEMENTS. 



B. Lattice models. 



A. Definition. 

We wish to use the concepts we have introduced in 
the preceding sections while allowing some motion of the 
grains, of small but finite extent, which might alter the 
list of closed intergranular contacts. Consequently, we 
introduce the assumption that displacements, from a ref- 
erence configuration, are small enough as to be regarded 
as infinitesimal quantities. This approximation of small 
displacements (ASD) is a crucial step that is very often 
taken in solid state mechanics. Indeed, it is indispens- 



Regular packings of monodisperse spheres in 3D (or 
discs in 2D) on FCC or hexagonal compact (respectively, 
triangular in 2D) lattices are simple systems that are of- 



truly monodisperse systems do not exist, and because 
of possible elastic deformations of the grains, one can- 
not expect such lattices to remain perfectly regular and 
undisturbed. However, as lattice perturbations will be 
small, it is a common practice @, Ufl [III [l?], [H, [29| to 
resort to the ASD, with a perfect lattice as the reference 
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configuration from which displacements and strains are 
evaluated. 

Consider e.g., the case of slightly polydisperse discs on 
a triangular lattice, as in systems A and B. A perfect 
lattice can be chosen as the reference state, in which the 
spacing between neighbouring sites is the lowest upper 
bound a of the diameter distribution. Diameters are as- 
sumed to be distributed between a(l — a) and a, with a 
small parameter a«l. The diameter of disc i is thus 

^ = 0(1-^), (29) 

Si being a random number, drawn independently for each 
i between and 1. When a certain number of intergranu- 
lar contacts is created, as it is often necessary (cf. section 
III) in order to sustain some external forces, the lattice 
will be slightly distorted, with displacements of order a. 
The ASD amounts to deal with all relevant quantities to 
leading order in a. In all possible contacts, the normal 
unit vector is kept parallel to one of the three directions 
of dense lines in the triangular lattice. It is convenient to 
work with a fixed structure So that comprises all bonds 
between nearest neighbours on the lattice. If grains are 
required to touch to exert a force on one another, forces, 
in a state of equilibrium under a supported load, will be 
carried by some contact structure, the bonds of which 
form a subset of Sq. 

One might then regard problem PA1, in system A, as 
defined on So- Once the random radii were fixed, we 
found, within the ASD, an equilibrium configuration for 
problem PA1, satisfying the Signorini condition 121} in 
which the contact structure was SA1. Similarly, once 
the values of the radii were known in system B, SB1 
was found, within the ASD, as the contact structure 
corresponding to a solution of problem PB1, posed on 
SB3=So. Within the ASD, all displacements and defor- 
mations are proportional to a, and the problem is, apart 
from a scale factor a for displacements, only sensitive to 
parameters (<5;)i<i<n- 

Such is not the case, of course, without the ASD, if 
one takes into account the rotations of unit vectors n; of 
the bonds due to the deformation of the lattice. 



V. ANALOGY WITH SCALAR PROBLEMS. 

We briefly recall the analogy between the mechanical 
problems we have been discussing, within the approxima- 
tion of small displacements, and that of current transport 
on a resistor network. Such an analogy was presented 
e.g., in ref. [20(. It is useful because some properties are 
more immediately intuitive in scalar models, and because 
statistical models (percolation, directed percolation, min- 
imum paths. . . ) have been more extensively studied and 
are more familiar in the scalar case. The term 'scalar' 
refers to the transport of a scalar quantity (current) as 
opposed to a vectorial one (force) in mechanical prob- 
lems. Currents entering one node by the conducting 
bonds of the network should balance the external current 



fed into that node, just like bond forces balance external 
efforts. The analog of the displacement vector (which, in 
the general case, also involves angular displacements) is 
the (scalar) potential of a node, and the duality between 
forces and displacements translates into the duality be- 
tween currents and potentials. All the developments of 
section II, adapted within the ASD to displacements in- 
stead of velocities, are valid for resistor networks. Sui is 
the potential drop in bond I. One may define spaces T , C, 
V, Sq, M operators G and G T , state the theorem of vir- 
tual power, etc. . . The analog of a system of self-balanced 
bond forces is a set of currents satisfying the conserva- 
tion law without any external source, i.e., a combination 
of current loops. One may define as many linearly inde- 
pendent elements of M as there are disconnected parts 
in the network. The number of degrees of freedom Nf is 
now equal to the number of nodes. It is related to the 
number of bonds N, the number of independent loops 
h and the number of disconnected parts (1 for a connex 
network) k by the scalar version of eqn. 1161 

N + k = N f +h, 

a simple topological identity valid for an arbitrary graph. 



VI. THE ISOSTATICITY PROPERTY. 
A. Statement and context. 

We consider an assembly of rigid, frictionless grains 
that only exert normal contact forces on one another. 
Those forces might however, be attractive or repulsive. 
We assume that the system, submitted to a prescribed 
load, has evolved to an equilibrium configuration in which 
the contact structure supports the load. We also re- 
gard the geometric definition of particles as incompletely 
known, thereby introducing randomness: such parame- 
ters as grain diameters or radii of curvature are to be 
regarded as distributed over small intervals. 

Then one can state the following remarkable property: 
with probability one, the problem, posed on the contact 
structure, is isostatic. 

Such an isostaticity property was (more or less explic- 
itly) reported in ref. [201 ] and articles cited therein, in the 
case of triangular lattice systems, within the ASD, with 
grains satisfying the Signorini condition |2"T1 Isostaticity 
was also stated in refs. 19. Il5l. fl7L [30j. that deal with the 
same model. Moukarzel [33l [34J then argued that sys- 
tems of frictionless grains interacting by repulsive elas- 
tic contact forces should become isostatic in the limit 
of large contact stiffnesses. And ultimately, Tkachenko 
and Witten [35[ derived an isostaticity property for disor- 
dered systems of rigid frictionless spheres in arbitrary di- 
mension, each grain being submitted to an external force 
(e.g., to its weight), whatever the sign of contact forces. 

Here, we will establish the isostaticity of the problem 
(h = 0), rather than the isostaticity of the structure 
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(h = and k — kg), in quite general situations. As 
we shall see in section VIII, full rigidity (k — /c ) in addi- 
tion to absence of hyperstaticity [h = 0), is a less general 
property, of geometric, as opposed to topological, origin. 



B. General arguments. 

The arguments we give below to establish the iso- 
staticity property emphasize the peculiarity of equilib- 
rium states, in which sufficiently many intergranular con- 
tacts should be created in order to resist the externally 
imposed forces. Thus such states belong to a subset of 
configuration space of vanishing measure. Grains have 
been brought to rest by some unspecified dynamic dis- 
sipative process. Our derivation admittedly retains a 
heuristic flavor, for a definitive proof would require much 
more specific mathematical assumptions. Readers that 
demand more mathematical rigour will have realized that 
arguments presented by other authors [33l 134 , |35( are not 
without reproach either, and may refer to the next para- 
graph. There, within the ASD (and thus at the expense 
of additional assumptions about the magnitude of dis- 
placements from a reference configuration), isostaticity 
is rigourously deduced. 

To ease the presentation of our arguments, let us intro- 
duce a few compact notations. We denote as (qi)i<i<N f 
a set of coordinates in configuration space £. The geom- 
etry of the grains depends on some random parameters 
(sizes, shapes. . . ), collectively denoted as £. £ might be 
regarded as a vector with a large number, say p, of com- 
ponents: C G The evolution of the granular system 
can be modelled as a function $ that maps an initial 
configuration ('?i)i<i<Ar / to the actual equilibrium con- 
figuration (qi)i<i<N f ■ The motion of the grains from 

(liW-liKNf t0 {li)i<i<Nf might e.g., be described by a 
differential equation. $ then expresses the dependence on 
initial conditions. <E> also depends on £, which has the role 
of a set of parameters. To proceed, on has to assume that 
this dependence is sufficiently regular: $ : £ x ]R P — > £ 
is generally a smooth function. Although the evolution 
of a pack of grains is expected to exhibit a high sensitiv- 
ity to parameters and initial conditions, it is dissipative 
and will bring the system very close to equilibrium in a 
finite time. Chaotic trajectories deviate fast from one an- 
other, but the evolution in a finite time is expected to be 
expressed by a smooth mapping, that also depends con- 
tinuously on parameters £, except perhaps for peculiar 
values that correspond to bifurcations between different 
sets of final states or 'attraction basins'. If, for instance, 
one reproduces the same dynamical evolution from the 
initial to the final configurations and gradually change 
the size of one particle, one expects, physically, the final 
state to change only gradually, until for some value of the 
geometrical change some rearrangement of finite extent 
will suddenly take place. We assume such bifurcations 
only occur for isolated values of the parameters, such 



that around the actual ( £ IR P , there exists generically a 
neighbourhood Q within which the parameter set might 
vary without creating any discontinuity or closing any 
additional contact in final configuration (?i)i<i<jv, € £ ■ 

Consider now the set L of intergranular contacts cor- 
responding to this configuration (the contact structure, 
as defined in section II). As £ changes within Q, main- 
tained contacts form some non-empty subset of L, which 
is sufficient to carry the load. 

If C £ ^ varies along a curve parametrized by u, so 
does, via the mapping <!>, (qi)i<i<N f in £ . If a contact 
(i,j) G L is to be maintained in this motion, one must 
have: 



This means that the coordinates of grains i and j have to 
adjust to the change in grain geometry £. If parameter 
u is formally regarded as time, relative normal velocities 

SVij — — , in all contacts that are maintained, are re- 

du 

quired to balance the effect of the change of to ensure 
that equality [30] is still satisfied. Increasing, if needed, 
the number p of £ components, it is natural to assume 
that such conditions on relative velocities are indepen- 
dent from contact to contact, for the required value of 
SVij only depends on those geometric parameters that 
govern the shape of grains i and j in the immediate 
vicinity of their contact point. Therefore, for a list L 
of N contacts to be maintained for arbitrary £ S f2, any 
A-vector (6Vi)i<i<n £ IR of possible relative normal 
velocities in the contacts of L must be compatible. In 
view of condition 1 1 51 only such contact structures L that 
are devoid of self-balanced sets of internal forces (i.e., 
such that h = or <So = {0}) can be maintained. If, 
exceptionnally, the equilibrium configuration (qi)i<i<N f 
admits one non-vanishing element (7z)i<i<jv of So, then, 
as the condition 

1<1<N 

cannot be ensured for arbitrary (SVi)i<i<n e M N , and 
grains cannot interpenetrate, one at least of the contacts 
I such that 7; ^ will open (SVi < 0) upon slightly 
tampering with geometric parameters £. 

We have thus shown that, with probability one, the 
contact structure in the equilibrium configuration cannot 
be hyperstatic, the degree of indeterminacy of forces h is 
equal to zero. 

The above derivation relies on rather specific assump- 
tions about mapping $. One should be aware, however, 
that we are free to choose any initial configuration that 
does not violate impenetrability conditions. The assump- 
tions we have relied upon are quite natural when the ini- 
tial and final equilibrium configurations are close to each 
other. Basically, one has then to accept the idea that the 
fine geometrical details of grain surfaces, in the vicinity 
of their contact points at equilibrium, do not significantly 
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influence their trajectories except in the very final stage. 
Thus they can be regarded as randomly chosen during 
this ultimate stage of the approach to equilibrium, as 
though the system 'realized' then what their actual val- 
ues are. In the next subsection it is assumed that the 
'initial' and final state are so close that the motion be- 
tween them might correctly be described within the ASD. 
Other derivations might resort to fictitious construction 
processes of the granular assembly, in which $ is replaced 
by a simpler function. One might consider, e.g., sequen- 
tially bringing the grains, one by one, to their equilibrium 
position, thus gradually enlarging the list of contacts. If, 
at any stage in the process, h is strictly positive, some 
of the contacts cannot be maintained on slightly altering 
some of the geometrical details of grain surfaces near the 
most recently created contacts. 

The equilibrium state, as we have just concluded, is de- 
void of hyperstaticity (h = 0). What about its possible 
mechanisms ? We have assumed that it can support the 
load. It is tempting to conclude that mechanisms do not 
exist in the generic case, since the orthogonality condi- 
tion [17] would have to be maintained as the shape of the 
grains is altered. However, one has to keep in mind that 
equilibrium configurations are very peculiar ones, and we 
shall see that the existence of mechanisms in the equilib- 
rium state depends in general on the sign of intergranular 
forces, and on the shape of the grains. 

C. Alternative derivation within the ASD. 
The special case of lattice models. 

A slightly different point of view may be adopted in 
the framework of the ASD: within the approximation, 
the problem being replaced by a simplified one, the iso- 
staticity property can be established in a rigourous way. 
Also, the analogy with the scalar problem might make 
the result more immediately intuitive. Let us assume the 
ASD to be valid with a reference configuration in which 
all contacts are slightly open: a list of bonds is defined, 
with strictly positive values of interstitial thicknesses h®. 
h®j , the distance separating the surfaces of grains i and 
j is to be regarded as a random number that depends on 
fine details of their geometry. /A values for the differ- 
ent bonds are independent and continuously distributed. 
Once the system has been brought to an equilibrium con- 
figuration, forces are carried by contacts, i.e. bonds I for 
which hi = 0. If (7/)i<;<tv is a set of self-balanced forces 
carried by those contacts, the theorem of virtual work, 
applied with such bond forces on the one hand, and with 
the displacements from the reference to the equilibrium 
configurations on the other hand, yields : 

5>(fc°-Jn)=£Ti*° = 0. (31) 

i=i i=i 

Thus a certain linear combination of the random dis- 
tances hf has to be equal to zero. Coefficients (7j)i<;<tv 



are fixed once the reference configuration is known. 
Moreover, via an iterative dilution process, they can be 
chosen among a finite set, as we now show: assume a set 
of self-balanced forces (ji)i<i<n to exist, and define the 
set Bo of bonds I for which 7; 7^ 0. Then, as long as it 
is possible, proceed to successive 'dilutions' of this set, 
defining B\, B2, etc. . . requesting that there is one bond 
less in Bk+i than in Bk, but that it is still possible to 
find self-balanced forces localized on the bonds of the re- 
duced set. The final Bk , that can non longer be diluted, 
will be such that the values of 7; will be uniquely deter- 
mined for each I € Bk , up to a common factor, which 
is fixed if one imposes the condition that the largest 7/ 
is equal to one. In this way, one thus defines irreducible 
sets of self-balanced forces, that are put in one-to-one 
correspondence with certain substructures of the whole 
contact structure. In a finite system, one thus has a finite 
number of such irreducible sets of bond forces. If a sys- 
tem of self-balanced forces can be carried by the contacts 
that are closed, then equation I5T1 has to be satisfied with 
one of the irreducible systems of self-balanced forces, an 
occurrence of probability zero. 

The scalar analog of this derivation is especially 
straightforward. To the requirement that only particles 
in contact exert a force on one another corresponds the 
condition that a bond between sites a and b on the re- 
sistor network can only carry a current when the poten- 
tial difference v a — Vb is equal to a prescribed value, v® b . 
Parameters v® b are to be regarded as random, chosen 
according to a continuous probability distribution and 
independent from bond to bond. Then, the appearance, 
once some current is injected at one node of the resistor 
network and extracted at another, of a loop of current- 
carrying bonds is to be discarded as an occurrence of zero 
probability. (One may of course define irreducible loops, 
as the ones that carry a unit current and do not contain 
stricly smaller subloops). Assume three bonds, making 
a loop between three sites, say 1 — > 2 — > 3 — > 1, to carry 
a non- vanishing current (figure [5]). This implies an exact 




FIG. 8: Three bonds forming a loop in the resistor network. 
A current might circulate as indicated by the arrows. 

relation of the form ±v® 2 ± «5 3 i "3 1 = 0> which has no 
chance to be satisfied. 

Let us consider now, as an example, returning to gran- 
ular systems, the small hyperstatic structure of fig. [51 and 
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assume the 7 grains have been brought, from the refer- 
ence configuration of the triangular lattice model defined 
in section IVB, in which all interstices are open (h®j > 0), 
to an equilibrium configuration in which the 12 bonds are 
closed contacts, with h^ — 0. Labelling the grains as on 
the figure, equation 1311 reads : 

7 6 

^ h>1,i ~ E ^M+l: 
i=2 i=2 

which is true with probability zero for continuously dis- 
tributed independent random numbers h®j. Within the 
lattice model with random diameters, as introduced in 
section IVB, one has 

h°ij = ^+6^, (32) 

one obtains a relationship between 5j's: 

Si = -(S 2 + S 3 + 5 A + 5 5 + S 6 + S 7 ), 
6 

which, once again, is satisfied with probability zero. 

It is less obvious, however, that the disorder on the 
radii of discs that remain exactly circular (or of perfect 
spheres in 3D) is sufficient, because of the induced dis- 
order on hij's, as in eqn. [32l to forbid the existence of 
any set of self-balanced contact forces. The problem is 
that, because of [32l interstitial thicknesses are no longer 
independent. On transforming [3"T1 into a relation between 
S^s, one gets 

e(e^U=°> 

which might well be satisfied if 7tt = f° r eacn i. 
This latter condition has no chance to be obeyed in a dis- 
ordered system, but may be achieved on a regular lattice. 
This does not occur, however, with nearly monodisperse 
discs on the regular triangular lattice in 2D, because 3 
independent conditions per disc are to be satisfied, and 
the number of contacts, at most three times the number 
of discs on this 6-coordinated lattice, has to be strictly 
smaller, because hyperstatic configurations like that of 
fig. [5] cannot exist. 

The situation is different for the analogous 3D model, 
defined with slightly polydisperse spheres on the sites of 
an FCC lattice. Each sphere has 12 nearest neighbours, 
and one may find hyperstatic structures in which con- 
tacts will be maintained with polydisperse spheres. A 
simple example of such a structure can be found, with 24 
spheres and 64 contacts 60] . Although a small amount 
of polydispersity eliminates hyperstaticity in 2D triangu- 
lar lattices of discs, it does not do so in FCC lattices of 
spheres, provided the grains, in spite of the distribution 
of radii, remain perfectly spherical. If the shape of the 
grains is also affected by the slight geometric disorder, 
then (with the notations of fig. HJ), one has ||Ry|| ^ ||R-ife|| 



for j =/= fc, interstitial thicknesses hij become independent 
in all bonds of the lattice, and hyperstaticity is forbidden. 
(Within the ASD, it is consistent to ignore the rotation of 
unit vectors due to small departures from sphericity) . 

D. Consequences. Remarks. 

Once the list of active contacts in an equilibrium state 
is known, isostaticity of the problem enables a purely geo- 
metric determination of the forces, independently of ma- 
terial properties. As an example, system C was brought 
in equilibrium under the load defined by eqn. 1191 with 
conditions [2TJ As soon as the list of contacts (structure 
SC) is known, the set of normal contact forces is entirely 
determined. 

This gives a meaning to the limit of rigid particles: 
in generic situations, when the sizes and shapes of the 
grains are affected by some amount of randomness, there 
is no problem of force indeterminacy once an equilibrium 
configuration has been reached. The actual value of con- 
tact forces will not depend on the detail of the contact 
law, provided it might be regarded as rigid, but it will 
be sensitive to fine geometrical details. As an example, 
consider frictionless elastic contacts obeying eqn.(2H Let 
us assume a stable equilibrium state of the grain assem- 
bly, regarding the grains as perfectly rigid ( condition l2"Tj) . 
has been reached. One thus has a local minimum of W 
(defined in eqns. [Hor^]). Then, let us take into account 
the finite, but small, deformability of the contacts. The 
same list of contacts will carry forces that, to first order 
in the small displacements, do not change. Evaluation, 
within the ASD, of relative normal displacements ft.; < 
in force carrying contacts yields hi = —(-^-) 1 ^ rn , such 
relative displacements are compatible because of the iso- 
staticity property, and the resulting elastic energy, 

w el = _J_y R \ h = _J_y R -l/m (m+l)/r 

TO + 1 ^ ' ' TO+1 ^ 1 H 

tends to zero as stiffness constants Ki tend to infinity. 
Thus the actual values of constants Ki and exponent m 
(these data might vary from contact to contact) are ir- 
relevant. 

Once an equilibrium state has been reached, force val- 
ues do not depend on the details of the contact law: this 
is an important step on the way to the reduction of the 
mechanics of granular systems to geometry-the basic goal 
of the present paper. This contributes to ease the deriva- 
tion of generic mechanical properties of granular systems. 

The simplification that results from the isostaticity 
property should however be balanced with the two fol- 
lowing difficulties. 

Firstly, configurations of granular systems, due to the 
same isostaticity property, are necessarily quite sensitive 
to fine geometric details: tiny variations of grain dimen- 
sions or positions might lead to opening of some contacts. 
As all contacts are indispensable to support the load, the 
system has to rearrange somehow to create other contacts 
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that compensate for one that were lost. This is the origin 
of a property known as fragility, to be more accurately 
defined, and discussed, in part IX. 

Secondly, one should be aware that the choice of an 
equilibrium configuration among several possible ones 
might depend on other physical parameters than the ge- 
ometry of the grains. The reduction to geometry is thus 
not complete. In section VII below, the consequences of 
the ASD are studied, and it is shown that mechanical 
problems are entirely geometric within the approxima- 
tion. 

As a consequence of the absence of hyperstaticity (h — 
0), one readily obtains, from ll6[ a bound on the number 
of contacts N that carry a force, involving the number 
Nf of degrees of freedom of the particles belonging to the 
backbone of the force-carrying structure: N < Nf — kg. 
Neglecting the effect of boundary conditions on the count 
of Nf in large granular systems, one gets an upper bound 
on the coordination number c = — : 

n 

J c < 2d for spheres ,„„n 

\ c< d(d +1) in the general case, ^ ' 

Particles in 3D that possess an axis of revolution, like 
spheroids, also have one trivial rotational free motion (in 
the absence of friction). Thus one should subtract one 
degree of freedom for each, hence the bound c < 10, 
instead of the general 3D value 12. 

Interestingly, an estimate c ~ 11 for the coordination 
number of long rods or fibers was given by Philipse (36| , 
on the basis of some statistical assumptions about the 
random packings of such particles. 

What we have established is in fact the absence of hy- 
perstaticity of a generically disordered assembly of rigid 
grains, regarded as frictionless. Forces, in the derivation, 
only appear as convenient auxiliary quantities ('virtual' 
forces) to deal with a purely geometric problem. The 
conclusions thus holds in the presence of solid friction. 
Assemblies of rigid grains with friction therefore abide 
by inequality 1331 (It is of course well known, from nu- 
merical simulations in particular [3^ . l38l . [3{|, that the 
contact coordination number is a decreasing function of 
the friction coefficient). 

It is also worth pointing out that [33] does not depend 
on the polydispersity of the grains. Grains that are much 
larger than their neighbours will often touch a large num- 
ber of them. However, this effect should be compensated 
in the average coordination number by an opposite one, 
affecting small grains. When they touch a large one, this 
latter effectively occupies half of the surrounding space, 
thereby reducing the possibility for other contacts. 

On the ground that force-carrying structures should 
be rigid (devoid of mechanisms, k = kg) the opposite 
inequality, N > Nf, whence the lower bound d(d + 1) 
(2d for spheres or discs) for the coordination number, is 
sometimes quoted in the literature [H, H(| • We regard 
it as wrong in general (although true for systems of non- 
cohesive rigid frictionless spheres, as we shall see). As 
pointed out by Alexander [40| . the physically relevant 



concept is not rigidity, but stability (under a given exter- 
nal load). This is discussed in section VIII below. First, 
section VII is devoted to the exploitation of potential 
minimization properties within the ASD. 

VII. EQUILIBRIUM AND POTENTIAL 
MINIMIZATION WITHIN THE ASD. 

The approximation of small displacements introduced 
in section IV has several important consequences. Find- 
ing an equilibrium state amounts, in some cases, to solv- 
ing a convex minimization problem, for which optimiza- 
tion theory provides useful properties and tools. The 
relationship with percolation or minimum path models 
are also to be discussed within the ASD. 



A. Convexity. 

When the potential energy is a convex function of dis- 
placements or positions, and when the rigid constraints 
define a convex set in configuration space, then the search 
for a stable equilibrium state is a convex optimization 
problem, and the following important properties can be 
exploited [HI ]. 

1. The equilibrium conditions, which express the sta- 
tionarity of the potential, are not only necessary 
conditions for potential minimization (i.e., stabil- 
ity), they are also sufficient. 

2. A local minimum of potential IT is a global mini- 
mum. W is flat, equal to its minimum value, over 
a convex set of possible equilibrium configurations. 

3. A structure being given, a supportable load will be 
supported. 

4. Equilibrium forces are the solution to another op- 
timization problem (the so-called dual problem). 

5. Rigid laws and elastic ones can be dealt with in the 
same way. 

Let us, among the contact laws presented in section 
III, distinguish the ones that lead to convex problems. 
It should be remarked first that standard convexity is 
defined in vector spaces, not on manifolds. In order to 
exploit the classical results of convex optimization the- 
ory to grains of arbitrary shape, it is necessary to place 
ourselves within the frame of the ASD, which replaces 
the curved configuration space by its flat tangent space 

As intergranular distances hi are, within the ASD, 
affine functions of displacements, it follows that both 
rigid constraints hi > or hi < define a convex set 
(and so does hi =0): the accessible part of configuration 
space is a simplex, a convex set whose boundaries are a 
collection of flat sections (parts of affine spaces). Since 
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the potential energy of external forces, W, is linear in 
the displacements, its minimization belongs to the class 
of linear optimization problems, that are the subject of a 
large literature in applied mathematics and operational 
research. This important case -granular systems within 
the ASD with contact laws of type ETJ or systems abid- 
ing by [23] or [24] , or tensegrities-is dealt with in detail in 
section VIIB. 

Still within the ASD, contact laws involving smooth 
interaction potentials will lead to convex problems if the 
potential function w is convex. This is the case for uni- 
lateral elasticity, as defined in [26] and [27] but not for in- 
tergranular potentials that possess an attractive tail like 
on figure [7] 

Outside the ASD, convexity can be discussed in the 
case of spheres or discs, since, ignoring rotations, their 
configuration space is flat. One immediately checks, 
then, that impenetrability constraints hi > 0, once hi is 
no longer approximated as an affine function of displace- 
ments, define a non-convex set of admissible configura- 
tions. The opposite inequality hi < 0, on the contrary, 
does lead to convex problems. As we shall see, friction- 
less spheres on the one hand, and systems of strings tied 
together on the other hand behave exactly in the same 
way, upon reversing the sign of forces and deformations, 
within the ASD, but strongly differ without the ASD. 



B. Rigid, unilateral contact law. 

1. Context. Notations 

The properties of convex problems enumerated above 
are valid, in particular, in the case of linear optimiza- 
tion problems, for which they are sometimes presented 
in particular forms [ill. |42| . Here, in order to stress their 
physical meaning, we shall directly rederive them. We 
consider an assembly of rigid frictionless grains, satis- 
fying the Signorini conditions 12 1 1 dealt with within the 
ASD. We assume a structure has been defined, and if the 
load is supported, some of its N bonds will, at equilib- 
rium, close (hi = 0) and transmit a force (fi > 0). The 
following also applies if condition [21] is replaced by l23l 
or M 

Keeping the same notations as in sections II and IV, we 
know that the impenetrability constraints are expressed 
with matrix G 

N f 

Forl<Z<A, ^GinUpKrf, (34) 

the transpose of which appears in the equilibrium equa- 
tions 

N 

Forl</x<iV/, Y. G ^=n Xt - ( 35 ) 
i=i 



Throughout this section, compact notations will be used 
for vectors of external forces (F ext for (F^ xt )i< fl <N f ) con- 
tact forces (f for (fi)i<i<N), interstices (hfor (ft;)i<z<jv), 
and displacements (u for (u a ,)k m <jv / ), the bracket nota- 
tion (e.g., (f |h)) is used for scalar products, while opera- 
tor notations and abbreviation for inequalities reduce [34] 
to Gu < h°. 



2. Minimization in displacement space. 

We now show that finding equilibrium displacements 
is equivalent to solving the following linear optimization 
problem: 

f Minimize W = -QX = -£ F/f 
1 1 with constraints: (|34|) 

We know from section III that a solution to problem V\ 
provides a set of Lagrange parameters (fi)i<i<N that sat- 
isfy both conditions [3l] and [35] (or |2"2"|) , and are therefore 
equilibrium forces. 

Conversely, in the case of a linear optimization prob- 
lem such as V\ , the stationarity condition is sufficient to 
ensure that W is minimized. 

This can be checked as follows: let u* S V represent 
one solution for displacements, and, likewise, let us de- 
note equilibrium contact forces as f * € IR^ . To u* corre- 
sponds the set of values h* for interstitial distances, and 
the Signorini condition might be expressed as 

(f*|h*)=0, 

while any displacement vector u € V, corresponding to 
h, satisfies 

(f*|h) > 0. 

From the theorem of virtual work, one then has 

W(u) - W(u*) = -(P|h*) + (f*|h) > 

and displacement u* minimizes the potential energy. 

Figure [5] is a schematic representation of problem V\ . 
A simplex, defined by a set of affine constraints like [34] 
is limited by flat faces, where some of the constraints 
are active. Its extreme points (the 'corners') are where 
a maximum list of constraints are simultaneously active. 
The criterion to be minimized is itself an affine function, 
it is constant on hyperplanes that are orthogonal to the 
load. Equilibrium is achieved on the simplex boundary, 
at least in one extreme point, in general on a simplex 
A in a space that is orthogonal to the load direction. 
Let k (smaller than Nf) denote the dimension of this 
space. Within the set of solutions, W is constant, and a 
certain number N* of contacts are maintained closed. Let 
us denote this structure as S*: it is the list of contacts 
that are closed for all equilibrium configurations. For 
those equilibrium states that are on the boundary of A, 
some additional contacts are created. It follows from its 



18 




FIG. 9: Aspect of simplex of variables satisfying affine con- 
straints like 1341 cut by the plane of the figure. W is constant 
on parallel hyperplanes (sketched as dotted lines, orthogonal 
to F, projection of load direction onto the plane) . W reaches 
its minimum at one extreme point at least (like A and B) or on 
'faces' or' edges', included in an affine space of dimension fe, 
that are part of the simplex boundary, (like segment AB). The 
hatched region is forbidden by impenetrability constraints. 



definition that k is the degree of velocity (here, within 
the ASD, of displacement) indeterminacy of S*. Since, 
from part VI, its degree of hyperstaticity is zero, one has 
k = N f -N*. 



3. Supportable loads will be supported. 

In general, displacements are thus determined up to 
some motion within convex set A. 

Let us now show that A is not empty if the load is sup- 
portable. We assume some statically admissible forces 
{fi)i<i<N to be defined on the bonds of the complete 
structure that was defined a-priori. Then a finite lower 
bound for W on the whole simplex of admissible displace- 
ments can be obtained upon writing the variation of W 
from the reference configuration as 

AW = - £ f?6m > - ]T 

1<(<JV 1<(<JV 

W, thus, cannot decrease to — oo within the simplex, and 
has to reach a finite minimum somewhere on the bound- 
ary. Moreover, one can show that A is also bounded, 
except for marginally supportable loads. We say the 
load is not marginally supportable if there exists a small 
neighbourhood of (F^ xt )i<^<N f in force space T within 
which all loads are supportable. Let us now consider 



a situation in which A is not bounded. One can then 
find one direction along which displacements go to in- 
finity within A. Now let us assume the load is not 
marginally supportable. One can apply a small load in- 
crement (SFffiiK^Nf, such that (F^ xt +5F^ xt )i< t ,< Nf 
is still supportable, with (SF^ xt ) 1 <^ l <N f in the direction 
for which A is not bounded, which leads to a contradic- 
tion. Therefore the load has to be marginally supportable 
if A is not bounded. 



4- Dual problem in bond force space 

We now turn to the dual optimization problem, to 
which equilibrium contact forces are the solution, viz. 

f Maximize Z(t) = -(f|h°) = - £j hffi (m) 
2 1 with constraints: (|35|) and f > ' 

We know that equilibrium displacements (u*) and con- 
tact forces (f * > 0) respectively satisfy [34] and [35] and 
are such that 

(f*|(G.u* -h )) =0. (37) 

Thus, any possible set of non-negative bond forces f bal- 
ancing the load is such that 

(f|(G.u* - h )) < = (f*|(G.u* - h )) 

on the one hand, and 

(f|G.u*) = (G T f|u*) = {F ext \u*) 

on the other, which entails Z(f) < Z(f*): f* is a solution 
to problem V2 ■ 

Conversely, if one starts from problem V2 , and consider 
a solution f * , then it is possible to define an 7V/-vector u* 
of Lagrange parameters corresponding to constraints 1351 
and an TV- vector h of non-negative Lagrange parameters 
corresponding to constraints f > 0, such that 

- h° + Gu* + h = 0. (38) 

Moreover, hi vanishes whenever // > 0. This means that 
u* is actually a displacement vector abiding by[34j and 
equation [35] entails that the Signorini condition, in the 
form [37] is a l so satisfied. We know then that u* is a 
solution to V\ . 

Equilibrium displacements and contact forces thus co- 
incide with the respective solutions to V\ and V2, a 
pair of linear optimization problems in duality. We have 
shown that: 

• If u* is a solution to V\ , then it is possible to find 
a solution f* to V2, E3 being satisfied. 

• If f* is a solution to V2, then it is possible to find 
a solution u* to V\ , [23 being satisfied. 
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• If u* and f * respectively abide by the constraints of 
optimization problems V\ and V2, and if, in addi- 
tion, [37] (equivalent to the Signorini condition !?!]) is 
satisfied, then u* and f * are respectively solutions 
to V\ and V2 ■ 

• The optimum value of the criteria are equal in both 
problems: condition [37] ensures W(u*) = Z(f*). 



5. The uniqueness property. 

Within the affine space of bond forces satisfying 135] 
constraints /; > define a simplex, and, just like for 
Vi, the set of solutions to V2 is a convex part B of its 
boundary. Let h denote the dimension of the affine space 
spanned by B. Since B is the set of possible equilibrium 
forces, h is in fact the degree offeree indeterminacy of the 
problem. Generically, one has, from part VI, h = 0, and 
the only solution to problem V2 is an extreme point of the 
simplex of admissible forces. We have thus shown that 
in terms of forces, the solution is uniquely determined. 
This is a stronger conclusion that the sole isostaticity of 
the problem established in part VI: in general, contact 
forces are uniquely determined once the list of contacts 
is known. In the case of a system of rigid grains, with 
contact law [21] dealt with within the ASD, the list of 
force- carrying contacts itself (the list of bonds, among 
those that are defined a-priori in the reference configu- 
ration, for which neighbouring grains will actually touch 
and exert a force on each other) is uniquely determined. 
Forces are carried by contact structure S* , which was de- 
fined in connection with the discussion of the solutions 
to problem V\, and, if some mechanisms exist (k > 0), 
the other contacts that might be created will not carry 
any force. 

If the contact law is [21] if geometrical changes from a 
reference configuration are small enough for the ASD to 
be valid, if the load is supportable (but not marginally 
so), then the system will reach an equilibrium state, 
which apart from bounded displacements within convex 
set A (that do not change W) is totally independent of 
all dynamical properties of the system, and entirely de- 
termined by the sole geometry. 



6. Examples. 

Systems A and B introduced in part II, were treated 
within the lattice model defined in section IV. B, with the 
ASD, and condition[2T] Structure SA1, once the random 
numbers Si were known, was obtained as the uniquely 
determined list of force-carrying contacts at equilibrium 
under the load defined byHHJ Within the ASD, it is pos- 
sible to close 2 other contacts, e.g., those that belong to 
SA2. However, they will not transmit any force. Like- 
wise, for specific values of the Si's, SB1 was obtained as 
the list of force-carrying contacts in system B submitted 



to the load that is represented on figure [2] It is possible 
to close some other contacts (such as those that belong to 
SB2), but they cannot carry (within the ASD) any force. 
Uniquely determined force-carrying structures, depend- 
ing on the load, will possess a varying degree of displace- 
ment indeterminacy k. Once system B, in addition to the 
forces on the perimeter, was submitted to small (ran- 
domly oriented) external forces exerted on each grain, 
then isostatic structure SB2 was obtained. 

In ref. @, the triangular lattice model, as in section 
IV. B, was studied for isotropic loads. As an application 
of the global minimization property, it was shown, within 
the ASD (to first order in a) that the maximum pack- 
ing fraction of polydisperse discs is, in the limit of large 
systems, equal to 

^ max = 27I (1 " fca) ' (39) 

with k = 0.314 ± 0.003 in the case of a uniform distribu- 
tion of radii. 



7. Minimal structures. Analogies with other problems. 

As equilibrium contact forces are the coordinates of 
an extreme point of the simplex of problem V2 , a max- 
imum set of inequality contraints fi > are simulta- 
neously satisfied as equalities, /; = 0. This means that 
force-carrying structure S* is minimal with respect to the 
equilibrium requirement 1351 In section VI. C, we invoked 
an iterative dilution process to define irreducible sets of 
self-balanced forces. Likewise, one can define minimal 
structures, such as S*, as irreducible by further dilution, 
since it is impossible to require more bond forces to van- 
ish if the load is to be balanced. Any such irreducible 
structure S might carry a unique set of bond forces bal- 
ancing the load, it geometrically determines one solution 
to equations 1351 

Recalling we have defined a loading parameter Q, to 
which all external forces are proportional, there exists 
for each minimal force-carrying structure S a set of coef- 
ficients (P[)i<i<n, such that the forces carried by S that 
balance the load are 

fi = PfQ- (40) 

By definition, one has 

f pf ^ if l e S 
\ (3? = if I i S 

Among all minimal structures S with non-negative coef- 
ficients j3f , S* minimizes 

ies 

Let us now recall the analogy with a problem of current 
transport on a resistor network, as introduced in part V, 
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with the following constitutive law. To the requirement 
that contact forces are repulsive corresponds an orien- 
tation of the bonds, which behave as diodes rather than 
resistors. Bond a — > b between nodes a and b carries some 
current i a b > that is related to the potential difference 
v a — Vb by the analog of the Signorini condition: 



tab > if v a - v b = v° b 
i ab = if v a - v b < v° b 



(41) 



The bond becomes a supraconductor (the analog of a 
rigid contact) when the threshold potential difference v® b 
is reached, and it is an insulator if v a — Vb is smaller. 

It is customary to define a scalar analog of the me- 
chanical load by injecting some external current I in one 
node, that we denote as i, and extracting it from an- 
other one, that we denote as o. / is then the analog of 
the mechanical parameter Q. A minimal structure (i.e., 
one that cannot be further diluted), to carry the cur- 
rent, is a path from i to o. If its coefficients (3 cannot 
be negative, it is a directed path, on which the current 
flow respects the a-priori orientation of the bonds. On 
such a path S, all bonds I € S carry the total current 
J, hence (3f = 1 for all I S S. In the analogous scalar 
problem, the current is carried by the directed path S* 
that minimizes, among all directed paths S from i to o, 
the criterion 



5>?- 

les 



In the scalar problem, the criterion reduces to a sum of 
'costs' associated with the bonds of the network. 

The analogous problem to V2 in the scalar case is thus 
the well-known minimum directed path (or directed poly- 
mer) problem on a network [43j |. This analogy was intro- 
duced in [2(j, for problem V\, upon transforming the 
minimum path problem into the dual problem, which 
consists in maximizing the potential drop Vi — v , knowing 
that in each bond / vi cannot exceed the threshold value 

. The dual point of view adopted here-the analogy for 
problem Vt~ stresses the geometric origin of equilibrium 
forces, as coefficients characterizing the maximum local- 
isation of efforts onto structure S*. Contact forces in 
granular packings have often been studied in the recent 
literature [l(J 0, EH . It is interesting to be able to de- 
fine them as the solution to a well-defined optimization 
problem of random geometry [ItJ ■ 

Some statistical properties of structures S* were stud- 
ied in refs. 0, [H| , in the case of the 2D triangular lattice 
model, as defined in section IV. B, with a uniform distri- 
bution of Si's. It was shown, in particular, for isotropic 
loads in the limit of large systems, that the density of 
force-carrying bonds tends to a non- vanishing limit, and 
the distribution of contact force values was evaluated. 

The statistical properties of the solution to the 'di- 
rected polymer' problem are related to those of directed 
percolation [43]. Likewise, one can expect, in the case, in 
particular, of a very wide distribution of values of ho in 



the mechanical problem, minimization problem Vi to be 
related to some unilateral percolation problem. Such a 
percolation model was never studied to our knowledge. It 
is a geometric problem, unlike generic central force per- 
colation [22| , for which (in 2D at least) only the topology 
of a diluted structure matters. 



8. Some macroscopic results for the triangular lattice 
model. 

To see what macroscopic mechanical behaviour might 
result from the properties stated in this section, we briefly 
recall here some results obtained by numerical simulation 
of the triangular lattice model [17| , as presented in sec- 
tion IV. B, with a uniform distribution of parameters Si 
(eqn.HU). 

Samples of up to 12600 discs were submitted to varying 
states of stress. The following inequalities, in which coor- 
dinate label 1 corresponds to one of the three directions 
of dense rows in the triangular lattice, and compressive 
stresses are conventionnaly positive, define the domain of 
supported loads, as macroscopically expressed in terms 
of stresses. 
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(42) 



All intensive quantities, like, e.g., distributions of force 
values, density of the contact structure, distribution of 
contact orientations, etc. ..were found to possess well- 
defined thermodynamic limits, independently of the de- 
tails of the boundary conditions, provided a uniform state 
of stress is imposed, and the stress tensor a satisfies con- 
ditions 321 as strict inequalities. Correlation lengths or, 
in other words, sizes of representative volume elements, 
or of independent subsystems, are finite, but appear to 
diverge as marginally supported loads (for which one of 
conditions H2l holds as an equality) are approached. 

Taking, as in section IV.3, the undisturbed lattice, in 
which the spacing between sites is equal to a, the maxi- 
mum disc diameter, as the reference state, a strain tensor 
e can be identified. It is related to displacement field u 
By 



du a 

dx/3 



dx a 



(43) 



and the potential energy per unit surface area is (sum- 
mation over repeated indices implied) 



W = 



(44) 



Coordinates of tensor e are found to be expressible as 
linear combination of the average of bond elongations Sui 
for the three bond orientations of the triangular lattice. 
In e space (3-dimensional for a 2D system), impenetra- 
bility conditions define, in the thermodynamic limit, a 
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strictly convex accessible domain V, limited by a smooth 
surface £, the equation of which we denote as 

fig) = 0, (45) 

while the interior of accessible region T> corresponds to 
the strict inequality: 

fig) < o. 

As a macroscopic consequence of the variational proper- 
ties stated in part VII, the relationship between tensors 
q_ and e is the following: 

df 

(Tij = A-A with A > , if /(e) = (46) 

aij = if /(|) < 

Wherever the granular system transmits stress, the value 
of e is as far as possible in the direction of a within T>, i.e., 
where the tangent plane to its boundary £ is orthogonal 
to (7, thus minimizing potential energy 1441 

T> is unbounded in the direction of non-supported 
loads. Strains go to infinity on surface £ when the stress 
tensor approaches one of the marginally supported di- 
rections. £ has three asymptotic planes, respectively or- 
thogonal to those three marginally supported load direc- 
tions. 

The one-to-one correspondence between supported 
stress directions on the one hand, and strain tensors such 
that /(e) = on the other hand, is a macroscopic trans- 
lation of the uniqueness property stated in paragraph 
VIIB5. The potential energy density has a finite thermo- 
dynamic limit (a result that generalizes to non-isotropic 
states of stress the one of equation I3U), and possible vari- 
ations of e within convex set A, discussed in VII. B. 3, 
shrink to a vanishing range (e becomes uniquely deter- 
mined) as the system size grows. 

Constitutive law |4"61 can be used to solve for stress and 
displacement fields whenever a sample of the model mate- 
rial is submitted to some external forces that do not lead 
to unbounded displacements and overall failure. The field 
of A values should be obtained on solving the full bound- 
ary value problem. 

C. Systems with bounded tensile forces. 

If the unilateral contact law [5U is replaced by (25J the 
remarkable properties stated above in VII. B are lost. Let 
us illustrate this on a simple example. Consider the sys- 
tem depicted on figure [101 to be dealt with, within the 
ASD, as a triangular lattice model in the sense of IV. B, 
the contact law being 1251 Only one disc is mobile (num- 
ber 1), and we first consider the case of a vertical force 
of intensity F y oriented downwards like on the figure, 
keeping F x = 0. (Later in part IX we come back to this 
simple example and discuss its behaviour when F x is al- 
tered). Two equilibrium positions are possible: disc 1 




FIG. 10: A small sample of the triangular lattice model, in 
which the only mobile disc, marked 1, is slightly smaller than 
discs 2, 3 and 4, and is submitted to an external force. Disc 1 
is shown in one possible equilibrium position, in contact with 
2 and 3. The other one is sketched with a dotted line. 



might either be in contact with discs 2 and 3, or with 3 
and 4. As grains are rigid and only exert normal forces 
on one another when they exactly touch, the problem 
is isostatic in both equilibrium configurations, in agree- 
ment with the general property of section VI. The load, 
defined with F y > 0, is always supportable on structure 
Si, consisting in bonds 1 — 2, 1 — 3, and it is also sup- 
portable on structure 5*2, consisting in bonds 1 — 3, 1 — 4 
as long as F y < foVS. 

Thus, for < F y < foV3, even within the ASD, the 
equilibrium state and the list of force- carrying contacts 
are not uniquely determined. Whether Si or S2 will be 
chosen depends on the trajectory of disc 1 from its initial 
(reference) position. 

Likewise, supportable loads are not necessarily sup- 
ported. To check this, let us remove disc 2. In its motion, 
disc 1 might come into contact with both 3 and 4, and, 
provided < F y < foV3, reach an equilibrium position, 
maintaining those two contacts. However, it might as 
well never meet disc 4, and find a trajectory, past disc 3, 
on which its potential energy will keep decreasing forever. 



D. Smooth, convex interaction potentials. 

In the case of the elastic contact law [26l within the 
ASD, all properties of convex problems enumerated in 
section VILA are valid. Let us state the 'elastic' ver- 
sions of the 'rigid' optimization problems of VII. B. Vi 
is simply replaced by 

Vf : Minimize W tot defined in|28J 



22 



while contact forces are the solution to 

vt | Maximize - £ [h^fl -^K^f^^ 
( with constraints: ([351) an d f > 

(47) 

The function of contact force / that appears within 
the sum is the opposite of the Legendre transform of the 
elastic energy w, regarded as a function of relative dis- 
placement Su, i.e., fSu — w(Su), taken with / = jj§^y- 

Thus, solving amounts to 'minimizing the comple- 
mentary energy', a common procedure to find the forces 
in an elastic problem. 

In fact, one could have defined a potential energy, 
in the rigid case, equal to +oo if grains interpenetrate, 
and treat rigid problems exactly like elastic ones, con- 
straint [34] being taken care of by the definition of the 
potential. If the region, in phase space, that is forbidden 
by the constraints is convex, then such a potential can 
still be regarded as a convex function. Both the condi- 
tion [2T] and elastic law [35] are then expressed by 

/ € dw(Su), 

in which dw(Su) denotes the subdifferential of w at Su, 
i.e., the set of all / such that w(Su') > w(Su)+ f (Su' —Su) 
for any Su' . This mathematical possibility to unify rigid 
and elastic laws is specific to convex problems. This is 
the precise meaning of property 5 cited in section VILA. 
Here, we preferred to resort to a separate presentation 
of the rigid case in section VII. B, to stress the physical 
consequences of the variational properties. The reader 
may refer to [46| for a more systematic approach. 

Comparing Vi and Vf , as defined by [35J and [17] one 
may expect the following behaviour for the distribution 
of contact forces, as a set of grains with elastic contacts 
is submitted to a constant load, but the stiffness con- 
stant K is gradually reduced. (Similarly, one could also 
increase Q, keeping K constant). When K is very large, 
the elastic term is negligible in comparison with Z(f ), and 
the values of the forces should coincide with the (unique) 
rigid contact solution of V% ■ Thus the contact structure 
should barely suffice to carry the load (isostatic prob- 
lem) , the forces should exhibit the characteristic disorder 
of granular systems, with large fluctuations, force chains, 
etc. . . On the other hand, let us assume that the list of 
possible contacts (structure So) is well-coordinated, that 
there are many more contacts that are easy to close upon 
increasing the confining forces or decreasing the contact 
stiffness parameters. Then, in the limit of small K, Z(f) 
will, in turn, become small in comparison with the elas- 
tic energy. The elastic term tends to share equally the 
forces between contacts. Thus, a narrow distribution of 
force values is expected in this limit, and spatial hetero- 
geneities should be strongly reduced. Knowing that the 
minimum structure S* and the complete list of possible 
contacts Sq are of comparable densities, the order of mag- 
nitude of the average force fo does not change as grains 



are made softer. The two extreme regimes of stiff and 
soft contacts should thus be respectively defined by the 

fo fo 
conditions K 3> — and K <C — , involving a typical 

K K 

interstitial distance ho. 

Those two limits, and the transition regime, in which 

the contact density increases, were observed jl6f on the 

2D triangular lattice model, as defined in section IV. B, 

with contact law |2"51 



E. Remarks. The 'elasticity' of rigid grains. 

As announced beforehand, we have exhibited, in this 
section, model granular systems for which, at the expense 
of several assumptions, including the validity of the ASD, 
mechanical properties are entirely determined by geom- 
etry. 

We have seen that the distinction between systems 
made of rigid or deformable grains is not necessarily as 
important as one might have expected: similar potential 
energy minimization properties might be stated, the limit 
of large contact stiffnesses might safely be taken without 
any singularity (subsection D), and macroscopic stress- 
strain relationships might be written for some systems of 
rigid grains, as recalled in paragraph B.8. The difference 
between the systems such that the search for an equilib- 
rium state is a convex minimization problem (in which 
case the properties listed in subsection A are satisfied) 
and the others, such as the example of subsection C, is 
finally more relevant. 

Constitutive law 25J expresses a one-to-one correspon- 
dence between the direction of stress tensor er and strain 
tensor e, which is restricted to belong to surface S. It is 
quite similar to a macroscopic elastic law, even though 
it applies to systems of rigid discs. The response to a 
supported stress increment will be reversible. If this in- 
crement, 5g_ is in the direction of the preexisting stress 
tensor er, then no additional displacement or stress will 
result for rigid grains. For deformable grains, if contact 
law [21] is replaced by 1261 a small deformation, inversely 
proportional to constant K, will follow. If, on the other 
hand, 5g_ is orthogonal to the initial stress tensor, its ap- 
plication will entail a small strain increment Se, such that 
the new strain tensor will be exactly the point of £ where 
the orthogonal direction is that of the new stress tensor. 
In this second case, the apparent elastic modulus is thus 
inversely proportional to the curvature of surface S. 

In spite of the analogy, presented in paragraph B.7, be- 
tween the backbone of the force-carrying structure and 
cost- minimizing directed paths for scalar transport, the 
statistical properties of those two systems are quite dif- 
ferent. In agreement with various results on disordered 
systems ofwains @, H3] the triangular lattice system was 
found d, [TBI, [ijj to possess a standard thermodynamic 
limit: intensive quantities like the density of the back- 
bone, the strains, the distribution of contact force values 
have limits in the limit of large system size (except for 
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marginally supported loads). On the other hand, unlike 
the force-carrying structure in the mechanical problem 
we have been studying, the optimal directed path in the 
corresponding scalar problem is a critical object. 

The validity of the ASD -that might at first sight ap- 
pear as a mere technical aspect- is finally a crucial ingre- 
dient of the model granular systems that we are studying 
here. The next section examines some stability proper- 
ties that are important as soon as one does not resort to 
the approximation. 



VIII. OUTSIDE THE ASD: QUESTIONS OF 
STABILITY. 

We now enforce, on physically acceptable equilibrium 
states, another requirement: that they should be stable. 
We limit ourselves to the cases when stability can be 
discussed in terms of a potential energy. If the equilib- 
rium state is a local minimum of the potential energy, 
then there exists a region of finite extent in displace- 
ment space, around equilibrium positions, within which 
the system is spontaneously attracted to the equilibrium 
configuration. 

Within the ASD, one can only discuss potential vari- 
ations that are of first order in displacements. When 
floppy modes exist (k > ko), they appear as marginally 
unstable and one cannot tell whether, to higher orders, 
they actually destabilize the equilibrium configuration. 
The mechanical response to small perturbations or load 
increments is strongly dependent on these stability ques- 
tions. 

In general, we will show, with examples (section A), 
that the answer might depend on quite specific geomet- 
rical features of the granular system, and on the contact 
law. We are only able to give general answers for spheres 
or discs, as shown in section B. Section C discusses some 
consequences on the geometry and coordination of gran- 
ular packings at equilibrium, and on the macroscopic me- 
chanical behaviour. 



A. Simple examples. 

We consider rigid frictionless particles of various 
shapes, and discuss the stability of simple configurations, 
that depends on the ability of contacts to withstand ten- 
sion, and on the shape of the grains. 

1. Bond alignments. 

Assume three spheres, or three discs in 2D, to have 
their centers aligned as on fig. [JTJ the two extreme ones 
being submitted to opposite forces in the direction of the 
line of centers. Let us discuss the problem in 2D. The 
determination of contact forces is an isostatic problem, 
and there is, apart from rigid body motions, a trivial 
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FIG. 11: An alignment of 3 spheres (left), the middle one 
touching the other 2. Spheres 1 and 3 are submitted to equal 
and opposite forces along the line of centers. The new equi- 
librium configuration, upon exerting a lateral force g on the 
midle sphere, is shown on the right. 

mechanism corresponding to free lateral motion of the 
middle disc 2. This is of course well known to lead to the 
familiar buckling instability if one pushes the extreme 
discs towards each other, and to be stable if one pulls on 
them, provided the contacts can resist tensile forces. In 
the latter case, assuming one controls the forces paral- 
lel to line 1-2 exerted on particles 1 and 3, while their 
position in the other direction is fixed, the system will 
respond elastically to a small additional force exerted on 
disc 2, even though the contact law is rigid. After the 
system reaches its new equilibrium state, the orientation 
of contacts is such that the new load is orthogonal to the 
floppy mode. Specifically, if g is the lateral force pulling 
disc 2 away from the line 1-3, and if / denotes the ex- 
ternal force exerted on 1 and 3, the new position of the 
center of disc 2 is such that, assuming equality of the 3 
radii, the angle 9o between 1-3 and 1-2 (fig. ITT]) is given 
by 

&o = tan- 1 (A), 
while contact forces (tensile, and therefore negative) are 

fl2 = f23 = -/COS(0 O ). 

The potential energy, as a function of 9 {9 parametrizes 
the free motion that maintains the two contacts), reads 

W = -2afcos(9)-agsm(9) = -a^Af 2 + g 2 cos(0-0 o ), 
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FIG. 12: An alignment like that of fig. 1111 the middle sphere 
being replaced by an object turning concave parts of its sur- 
face towards spheres 1 and 3. 

and has its minimum for 9 = 9q. 

This elastic behaviour is similar to that of a rigid string 
under tension, which will deform in response to lateral 
sollicitations. 

On carrying out the same calculations in the case of 
compressive forces, with / < 0, one will notice that g 
and 9q, corresponding to the equilibrium position of disc 
2, are now of opposite signs. One then has 

W = ay/Af 2 + g 2 cos(6> - 9 ), 

which is maximized in the unstable equilibrium position 
9 = 9 . 

In section VIIIB, we show that the conclusions reached 
on this simple example are general: any floppy mode in a 
system of discs or spheres that admits only compressive 
contact forces leads to an instability. If, on the contrary, 
all contact forces are in fact tensile, the system being thus 
analogous to a network of tight strings, any floppy mode 
is stable, and an elastic response to small load increments 
can be observed. 

Let us now replace disc 2 by a particle presenting con- 
cave surfaces toward discs 1 and 3, as shown on fig. [T2J 
The system is similar to that of fig. [Til the free lateral mo- 
tion of the middle particle, maintaining the contacts, is a 
mechanism. It is not difficult to show, however, that the 
configuration of fig. [12] has, compared to the alignment 
of discs, opposite stability properties: the mechanism is 
stable for compressive forces, unstable for tensile ones. 
Thus stability properties are quite sensitive to particle 
shape. 

2. Arches. 

Systems submitted to gravity provide other familiar 
examples of non-rigid equilibrium states. A string of cir- 
cular, or spherical, particles, each of them tied to two 
neighbours by a frictionless contact condition that sup- 
ports tension, behaves as a chain, and will eventually 



adopt a stable equilibrium configuration if one fixes its 
two extremities and let it dangle under its weight. The 
number of mechanisms in this system is equal to the num- 
ber of free particles. 

The analogous system to the chain, in which contacts 
transmit compressive forces, is the arch, fig. [T2J The 
general result for spheres entails that all arches made of 
spheres are unstable. However, one usually builds arches 
with appropriately shaped stones, e.g., carving them to 
share common flat lateral surfaces with their neighbours, 
as on fig.[T3T Such an arch is a system that possesses one 




FIG. 13: An arch built with stones sharing flat lateral sur- 
faces. 

floppy mode per stone (still assuming no friction), but 
its geometry might be adequately chosen to support the 
load. In such a case, any free motion of the stones, that 
slide on their flat common surfaces, all contacts being 
maintained, does not change the potential energy. One 
thus has an example of marginal stability. Such an arch is 
only able to carry the one particular load for which it was 
specifically designed. (Any amount of friction, however, 
stabilizes the system). 



3. A stable mechanism with strictly convex cohesionless 
grains. 

In view of the previous examples, one might be 
tempted to infer that mechanisms, when contacts only 
support compression, can be stable with concave grains 
(fig. I12p. are sometimes marginally stable with flat sur- 
faces (fig.[T3J|, but are always unstable with stricly convex 
grains (fig. Hip . This is however not true, as shown by 
the simple example of fig. [JJ] We are not aware of other 
general answers to this question of stability than the ones 
that are given for spheres below. 
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FIG. 14: The upper grain 1 relies on two of its neighbours and 
is submitted to its weight, oriented downwards. Its rotation 
is a stable floppy mode. 



B. General results for spheres and discs. 

1. Tensile contact forces (systems of cables). 

In the case when all contacts, at equilibrium, carry a 
tensile force, then stability is immediately proved once it 
is realized, as remarked in section VIIA, that minimizing 
the potential energy is a convex optimization problem 
(see property 1 stated in section VIIA). 

Just like for the simple example of figure [TTJ floppy 
modes can exist in stable equilibrium configurations. 
Then, the system will respond elastically to small load 
increments that provoke small motions of those floppy 
modes. Applying such load increments amounts to 
slightly deform the potential energy landscape on the 
manifold of configurations that maintain the initially ex- 
isting contacts. A new minimum is found, close to the 
previous one. 

Systems of rigid cables, whatever the level of defor- 
mation, should therefore possess exactly the same kind 
of elasticity, due to preexisting stresses, as assemblies 
of rigid frictionless particles without cohesion within the 
ASD (whose mechanical response to load increments was 
discussed in section VILE). 

Those properties were in fact discussed by Alexan- 
der [io| . in his monograph on the elasticity of various 
kinds of networks and amorphous systems, in the case 
when the contact law is elastic. Alexander pointed out 
that stable configurations are not necessarily rigid. He 
stressed that force-carrying bonds or contacts always 
have a stabilizing effect when they transmit a traction, 
and a destabilizing one when they transmit a compres- 
sion. 

Our present study, in this subsection, might be re- 



garded as complementary to his, since we deal with rigid 
contacts. 



2. Cohesionless grains. 

Let us now show that, in the absence of tensile force 
in the contacts, an equilibrium configuration of rigid, 
frictionless discs or spheres is necessarily unstable if the 
backbone is not a rigid structure. 

We shall do so by yet another application of the theo- 
rem of virtual power, as follows. 

We assume a packing of spheres to be in equilibrium 
under a prescribed load. Spheres are rigid, and the prob- 
lem is therefore isostatic, h — 0. Flat walls can also exist, 
e.g., as a device to enforce some kind of boundary condi- 
tion on the packing, but we assume that they cannot ro- 
tate. We assume there is at least one mechanism: k > 1. 
Consequently, it is possible to move the grains (and the 
walls) while maintaining the whole list of contacts. (The 
possibility that a mechanism could exist for the consid- 
ered equilibrium configuration alone, and disappear as 
soon as the grains are displaced is to be discarded as 
non-generic. This would, in particular, due to entail 
h > 1). We now study the variation of the potential en- 
ergy in one such motion, with a 'time' t parametrizing 
the trajectories, and show that it decreases. 

Objects do not rotate in this motion (this is an as- 
sumption for walls, and rotations of frictionless spheres 
are ignored anyway). Particle i has a time-dependent 
velocity Vj(i), and initially, in the equilibrium configu- 
ration from which the motion starts at t = 0, touches 
its neighbour j in a point A^., where the normal unit 
vector to its surface, pointing to the center of j, is , 
the equilibrium contact force being fy. Let Aij(t) de- 
note the material point of the surface of grain i that was 
at A®j initially. Similarly, following the material motion 
of j, one defines Aji(t), which does not coincide in gen- 
eral with Aij(t). It is possible, at each time t, to apply 
the theorem of virtual power, thus evaluating W'(t), the 
time derivative of potential energy W at time t, as fol- 
lows. The definition of a structure, in part II, was in fact 
completely arbitrary. Here, let us use this one: at time 
t, although objects i and j that are in contact effectively 
touch each other by a different point, define a bond to 
exist between Ay(i) and A^(t), oriented by n^, which, 
because objects do not rotate, is still carried by the com- 
mon normal direction to the surfaces of i and j at these 
two points. This structure might be used to define vir- 
tual, fictitious bond forces, that we choose equal to the 
initial equilibrium contact forces, i.e., fy, carried by n^- 
in the bond between Ay(t) and Aji(t). These forces are 
now used in the theorem of virtual work, with the real 
velocities. This is perfectly valid, because for each t 

• the virtual internal forces balance the constant load 

• in the bond between i and j, the force exerted on i 
is still equal to the opposite of the force exerted on 
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3- 

One obtains: 

the sum running over all bonds. As fijti^j does not de- 
pend on t, this is easily integrated. Denoting as Uy (t) 
the vector of origin Aij(t) and extremity Aji(t), the net 
variation of potential energy at time t, from the begin- 
ning of the motion is 

W(t) - W(0) = J2 /ynO.Uyft). (48) 

i<j 

In the motion, Aij(t) and Aji(t) are still extreme points 
of solids i and j in the respective directions n^- and — n?-. 
As spheres z and j have stayed in contact, it follows that, 
as shown on figure [TBI the contribution of bond i—j tol4"81 
is strictly negative, unless Aij (t) = Aji (t) , in which case 
it is zero. The same conclusion holds true for a contact 




FIG. 15: Sketch of the position, at time t of two spheres in 
contact. 

between a sphere and a flat wall that does not rotate. 
Consequently, one must have: 

W{t) - W(0) < 0, 

unless all intergranular contacts that carry non- vanishing 
equilibrium forces are maintained, in the motion, via the 
same material points. This latter condition means that 
the backbone of the contact structure in the equilibrium 
configuration moves as a rigid body. 



Mechanisms that only affects grains that do not carry 
any force, without altering the geometry of the backbone 
will not, of course, change the value of W and lead to 
instabilities. 

Otherwise, the instability is always present. We have 
shown that the backbone of the contact structure, in a sta- 
ble equilibrium configuration of a packing of rigid, fric- 
tionless spheres that do not support tensile forces in the 
contacts, is devoid of mechanisms other than rigid-body 
motions: k — ko- As we already knew, from part VI, that 
it cannot possess self-balanced contact forces (h = 0), one 
reaches the conclusion that it is an isostatic structure. 



C. Consequences. Discussion. 

1. Coordination of packings. 

The isostaticity of the force-carrying structure in pack- 
ings of rigid frictionless spheres with contact law !21l thus 
results from a stability analysis. The opposite inequality 
to the ones established in section VI. D, can, in this case, 
be stated: one has N > Nf, and consequently, N = Nf, 
on the backbone of the contact structure. For large sys- 
tems, the absence of floppy mode implies a lower bound 
on the coordination number: 

c> 2d on the backbone. 

This is equal to upper bound 1331 hence the equality: c = 
2d. 

However, for frictionless grains with different shapes, 
or for spheres with cohesion, one cannot expect in general 
inequality [33] to hold as an equality, even on the sole 
backbone. 

Returning to cohesionless packings of spheres, when 
each one is submitted to an external force, it has to be- 
long to the force-carrying backbone, and the whole sys- 
tem satisfies N — Nf (or, asymptotically for large sizes, 
c = 2d). This happens in system A, treated without re- 
sorting to the ASD. The force-carrying structure that was 
obtained, SA2, is isostatic and spans the whole system. 
When external forces are transmitted from the bound- 
ary, as in system C, floppy modes can exist, typically as 
isolated spheres, like discs 10 and 14 on fig. [3] or small 
sets of spheres, that are not or insufficiently connected 
to the backbone. If not too widely polydisperse sys- 
tems of spheres, regions that are totally shielded from 
force transmission are usually quite small. According 
to our experience in numerical simulations, if the radio 
of the largest to the smallest radius is 2 in a polydis- 
perse assembly of discs, then one very rarely sees more 
than 3 discs together in such regions. In 2D, ring-like ar- 
rangements surrounding discs that carry no force, such as 
29-30-31-15-6-5-13-28 and 11-3-9-22-23-24 
on figE] cannot easily be made very large: the curvature 
of the 'ring' would then decrease, increasing the risk of 
inward buckling. 
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2. Lattice models with and without the ASD. 

The triangular lattice model, as denned in section 
IV. B, of which systems A and B are particular samples, 
provides vivid examples of the difference between tensile 
contacts (systems of strings, satisfying |2"5|) and compres- 
sive ones (rigid grains obeying [21]), once dealt with out- 
side the ASD. Within the ASD, both types of systems 
share the same properties, and an equilibrium state of 
one of them can be mapped onto an equilibrium state of 
the other, as follows. In the reference state, rigid discs 
do not touch, since h% = (a/2)(#j + Sj)a > 0. This 
can be mapped onto a string network system, in which 
the 'contact law' is [221 on replacing each Si by — Si and 
attributing the length a(l + a(Si + Sj)/2) to the string 
joining i and j. On reversing the sign of external forces, 
an exact correspondence is achieved between equilibrium 
states. 

Fig. [TBI shows the force-carrying structure, as obtained 
within the ASD, in a hexagonal sample (for one random 
choice of Si values, drawn according to a uniform distri- 
bution) of 1141 discs. This system is submitted to an 
isotropic pressure, via a imposed homogeneous shrinking 
of the perimeter. 




FIG. 16: Triangular lattice model, within the ASD: force- 
carrying structure S* in a hexagonal sample submitted to an 
isotropically compacting load. Line widths are proportional 
to force intensities. The very same structure is observed in a 
corresponding system of strings undergoing isotropic tension. 



As established in section VIIB, such a structure is, 
within the ASD, only dependent on the random parame- 
ters Si . The dynamics ruling the motion of the particles 
from the reference to the equilibrium positions, and the 
actual value of a are both irrelevant. In the correspond- 
ing system of strings submitted to isotropic tension, ex- 



actly the same force pattern is obtained at equilibrium. 
We denote as S* the backbone of the contact structure, 
as displayed on fig. \TEl Just like in structure SB1, which 
carries the force in a similar sample of smaller size, many 
discs do not belong to S*, which only contains 619 of 
them, thus possessing 1239 degrees of freedom (counting 
the one of the 'wall'). Many floppy modes are present, 
381 of them are associated with bond alignments (discs 
having two contacts in opposite positions), and the re- 
maining 5 are more collective (like the one of figd]). Some 
statistical properties of S structures in the large system 
limit were studied in [lil ]. 

We numerically determined force-carrying structures 
in the rigid disc system under compression, and in the 
corresponding system of strings under tension, without 
the ASD. Those structures, that were obtained with a — 
1/48 (this value is now relevant), are respectively denoted 
as SC and ST, and shown on figures [T7l and [T51 Slight 




FIG. 17: Structure SC that replaces 5"* ouside the ASD in 
the case of contacts resisting compression. 

distortions of the regular triangular lattice, although not 
apparent on the figures, were taken into account in the 
calculations. From part VII, we know that ST is still 
determined by the sole system geometry: since forces 
are the solution to a convex optimization problem, the 
uniqueness property still holds. This is not the case for 
SC, and the result now depends on the actual dynamics 
(the rule that was adopted to move the discs to their 
final equilibrium positions). The calculation was carried 
out with the 'lubricated granular dynamics' method of 
refs. 0,[Hj]. 

As expected, SC is devoid of mechanisms: it is an iso- 
static structure, with 1052 discs, 2105 degrees of freedom, 
and exactly 2105 contacts. Only 89 grains out of the to- 
tal number 1141 do not belong to SC. Most of them are 
isolated grains, or pairs of neighbours (slightly larger re- 
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FIG. 18: Structure ST that replaces SC for an analogous 
system of cables resisting tension. 



gions shielded from the forces appear near the perimeter, 
due to a boundary effect). 

On the other hand, ST stays more tenuous, with 840 
discs only, and 1401 contacts. Thus 280 floppy modes 
still live on ST, 232 of which arc simple bond alignments 
and 48 are collective. 

In spite of those differences between the density of S* , 
SC and ST, it does appear on the figures that the spatial 
distribution of the forces is very similar, the strongest 
'force chains' remaining unaltered. The distributions of 
force values in S* and SC, in the limit of large systems 
were evaluated in ref. [l5| , and shown to coincide, within 
statistical uncertainties, except for the small forces that 
appear on SC in the additional contacts created by the 
buckling instabilities in S*. 

Thus, resorting to the ASD is quite a legitimate pro- 
cedure, provided a is small enough as to allow to regard 
the differences between SC or ST on the one hand, and 
S* on the other, as refinements that can be neglected. 

In the limit a — > 0, any contact force on SC is expected 
to tend to its value in S* , although the density of force- 
carrying contacts is discontinuous. 

In the system of strings under tension, on the other 
hand, mechanisms do not lead to instabilities, and the 
density of the backbone itself should continuously ap- 
proach that of S* as a — > 0. 



3. Role of grain shape: are spheres special ? 

We have seen that it is necessary to examine, beyond 
the ASD, questions of stability, to find qualitative differ- 
ences between intergranular contacts that resist compres- 
sion and cables that resist tension, and between spheres 



and other shapes. 

Of course, one expects macroscopic properties of gran- 
ular assemblies to smoothly depend on grain shape: pack- 
ings of nearly spherical grains will resemble packings 
of spheres. Experimentally, it has sometimes been ob- 
served that systems of spheres, in a quasi-static experi- 
ment, yield particularly noisy responses. It is also em- 
pirically known in civil engineering that granulates made 
of smooth and rotund particles, like river-bed gravel, are 
especially unstable and prone to large plastic deforma- 
tions. 

Unfortunately, detailed data at the microscopic level 
on non-spherical grains close to equilibrium are scarce. 

Although detailed analyses of such features are lack- 
ing, and our study of granulate stability should be ex- 
tended to the case of spheres with friction, one might 
speculate that such particular behaviours of rotund ob- 
jects could be related to the specific property we have 
established here: whenever some motion is smoothly ini- 
tiated (i.e., with a very small initial acceleration), while 
existing force-carrying contacts are maintained, then it 
will entail some loss of potential energy, and thus accel- 
erate further. Hence probably the jerky aspect of system 
trajectories in configuration space. 

Section IX discusses, precisely, when and how a system 
jumps from one equilibrium state to another. 



IX. MECHANICAL RESPONSE TO LOAD 
INCREMENTS: TOWARDS MACROSCOPIC 
BEHAVIOUR. 

So far, we have mainly dwelt on mechanical proper- 
ties of model granular systems. Those can be proved 
directly. We wish now to discuss possible macroscopic 
consequences in terms of the constitutive laws that are 
relied upon in a continuum mechanics description. Wc 
thus have to infer some of the properties of granular pack- 
ings in the limit of large systems. To be quantitative, 
some statistical knowledge of the geometry of large gran- 
ular systems is needed, which requires experiments or 
numerical simulations. Here, as we do not present new 
experimental or statistical studies, we shall focus on qual- 
itative properties, extrapolating on the characteristics of 
finite systems we have been presenting so far, and ex- 
ploitin g so me recent numerical results, especially those 
of ref. [r7j| . recalled in paragraph VIIB.8. 

Some macroscopic aspects of granular mechanics are 
recalled in part A. Possible origins of plasticity are dis- 
cussed in part B, in relation to grain-level characteristics. 
Part C examines some consequences of the strong iso- 
staticity property of systems of frictionlcss spheres with- 
out cohesion, in which case some response functions to 
load increments are related to the operator G, defined 
in section II in relation to equation 1121 corresponding 
to the isostatic structure. Part D exploits the results 
of ref. [ItJ , deriving the form of the macroscopic equa- 
tions to be solved when a small load increment is applied. 
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Finally, these results are compared, in part E, to some 
other approaches and theories, that were put forward by 
several authors in the recent literature, both at the mi- 
croscopic [H, H, UK and the continuum [H E El Q 
level. 



A. Macroscopic granular mechanics: known 
features, conflicting models. 

A classical way (see, e.g.,, in |l9() to study the macro- 
scopic mechanics of granulates is to submit a sample to 
a triaxial test. Such a device is designed to impose a 
uniform state of stress throughout the sample. It does 
not matter, for our discussion, whether this macroscopic 
stress is imposed via a fluid pressing on a flexible mem- 
brane (as in a laboratory apparatus, for the lateral con- 
finement) or via a control of the position of a rigid wall 
(as in some numerical simulations). We just need to re- 
member that a varying load is imposed, and depends on 
two parameters p and q, the axial stress {o~ yy on the fig- 
ure) being equal to p + q and the lateral one (a xx ), to p. 
A typical experiment consists in gradually increasing q at 
constant p. One may then observe the resulting strains. 
The classical elasto-plastic constitutive laws that are ap- 
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FIG. 19: The triaxial experiment. 

plied to granular materials are incremental, which means 
that they do not relate stresses and strains directly, but 
predict the increment of strain resulting from an incre- 
ment of stress, given the current state of the system (the 



definition of which might require other, 'internal' vari- 
ables) . Cycling sollicitations of small amplitude usually 
yield, in the stress-strain plane, loops with some amount 
of hysteresis. The surface area of such a loop as OABO 
on fig. [20] is the plastically dissipated energy associated 
to deviatoric stresses (to which the work due to volume 
changes has to be added to get the total plastic work). 

In marked contrast with classical soil mechanics ap- 
proaches, some authors recently proposed a new type 
of macroscopic mechanical description for the statics of 
granular packings [lj, [HI, 0, Ejj. According to them, 
resorting to strain variables should be avoided and one 
should look for direct relationships between the compo- 
nents of the stress tensor, so that it is possible to deter- 
mine the whole stress field in a granular sample by solv- 
ing hyperbolic second-order partial differential equations. 
Those, like wave equations, possess characteristics, pre- 
ferred directions along which they reduce to simpler, first 
order forms. To solve the problem, one may integrate 
along the characteristics that emerge from every point 
where some external force is applied. Consequently, in a 
packing in which the forces exerted on the top boundary 
(wall or set of particles) are known, a perturbation (exter- 
nal force increment), will propagate downwards, but will 
not be felt above the point where it is applied. The exact 
relation between stresses to be used should then depend 
on the actual process by which the sample was made. If 
the current stress level is changed, by, say, a manipulation 
of the boundary conditions, like in the triaxial test, then 
the granular system rearranges until the new constitutive 
relation, corresponding to its new state, agrees with the 
new externally imposed stress values. Those theories, in 
their current state of development, do not predict the ex- 
tent to which the system has to rearrange, or, in other 
words, the magnitude of the ensuing strain increment. 
It has been recently proposed [H[ that isostaticity could 
justify such theories for frictionless assemblies of grains. 
These suggestions are discussed in section IX. E below. 

We now turn to a discussion of some possible micro- 
scopic origins of plastic dissipation. 



B. Origins of plastic dissipation. 

When a given supported external load places the sys- 
tem in a uniquely determined equilibrium state, one has 
to expect a mechanical behaviour devoid of plastic dissi- 
pation. Hysteresis loops like those of fig. [20l cannot occur. 
Plasticity is related to the lack of uniqueness of equilib- 
rium states. At the level of continuum mechanics, it is 
sometimes termed 'internal friction', since the material 
behaves as if different layers of matter slided, with fric- 
tion, on one another within the bulk of the sample. We 
have thus identified two microscopic origins of internal 
friction in systems of frictionless grains. 

1. bounded tensile forces in the contacts (as in section 
VI.C) 
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FIG. 20: Schematic aspect of response to cyclic variations in 
q in the £22 - q plane. 
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FIG. 21: Loading parameter Q = versus coordinate x of 
the mobile disc of the system of fig 



to 2, 3, and 4, by cables that are slightly longer than the 
common distance between 2 and 3, and between 3 and 4. 



2. rearrangements of finite extent {i.e., the ASD is no 
longer valid) between equilibrium position of as- 
semblies of spherical grains. 

Let us illustrate these different behaviours on the sim- 
ple example of fig. [TO] (section VI. C). 

Starting from an equilibrium configuration in which 
the external force on disc 1, in contact with 2 and 3, is 
vertical, let us gradually increase its horizontal compo- 
nent F x . We first discuss the problem within the ASD. 
It is then a particular example of Vi discussed in section 
VII, a linear optimization problem with two unknowns 
(the coordinates of disc 1). In fact, the simplex within 
which potential energy W has to be minimized is exactly 
the one that was shown on fig. Points A and B on 
that figure are respectively the equilibrium positions of 
the center of disc 1 when it is in contact with 2 and 3, 
and with 3 and 4. Changes from one position to the other 
happen when the direction of F is orthogonal to that of 
segment AB. One may monitor the abscissa of the mo- 
bile disc, x, which, as presented on fig. [2T1 is related to 

loading parameter Q = -5 s - via a step-like function. 

y 

In analogy with this problem of rigid grains, one may 
build a system of rigid cables (resisting tension, but not 
compression), which, if treated within the ASD, yields 
exactly the same simplex of accessible configurations, the 
same optimization problem (Vi) as that of fig. [9] This 
system of cables is shown on fig. Node 1 is now tied 




FIG. 22: System of cables, equivalent, within the ASD, to 
the system of discs of fig. 1101 with the same values of external 
forces. Here, the cables joining 1 to 2 and 3 are taut, while 
the one joining 1 to 4 is not. 

Outside the ASD, the potential minimization problem 
for the system of cables is no longer a linear optimiza- 
tion problem, but, according to the general properties 
discussed in section VIII, is still a convex problem. In 
the plane of the coordinates of node 1, the simplex of 
fig. [9] changes into a domain limited by curved faces, as 
shown on fig. [23] The curvature of the faces being ori- 
ented inwards, this domain of accessible configuration is 
convex. When the orientation of force F is such that, on 
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FIG. 23: Minimization problem in the plane of coordinates of 
node 1, for the system of figure [22l without the ASD. The ac- 
cessible part of configuration space (outside the hatched zone) 
is convex. Its boundary has sharp corners (A and B), but, 
unlike on fig. [9] corresponding to the same problem within the 
ASD, displays curvature in between. Tangents to that curve 
at A and B are drawn. 



fig. [23] the direction of constant potential energy lines 
lies between those of tangents to the accessible domain 
in A and B, the equilibrium position is a point on arc 
AB, and only one cable is taut, the one joining 1 to 3. 
In this case, the motion along arc AB is a mechanism, 
but stability is maintained, just like in the example of 
fig. [TTJ There is still a one-to-one correspondence be- 
tween Q = |f- and x : as shown on fig. [24] As the differ- 
ence between cable lengths and distances 2-3 and 3-4 de- 
creases, displacements get smaller and smaller. The dif- 
ference Qb — Qa tends to zero, the curvature of the acces- 
sible region boundary on fig. 1231 vanishes, and the curve 
of fig. l24l approaches the ASD case, fig. [2T] Over a finite 
interval between Qa and Qb, the force-displacement re- 
lationship is a smooth function, unlike the stepwise de- 
pendency shown on fig. [21] (corresponding to the limit of 
very small motions). 

Let us now deal with the system of fig. [10] (with rigid, 
impenetrable discs and frictionless contacts that do not 
resist tension) outside the ASD. The accessible domain 
in the coordinate plane is, as opposed to the previous 
cases, no longer convex, as shown on fig. [25] The upper 
limit Qa of the Q interval for which position A is stable 
is now larger than the lower limit Qb of the Q interval 
for which position B is stable. Because of this bistability 
for Qb < Q < Qa, the Q versus x relation now exhibits 
hysteresis, as shown on fig. [26] 

As shown in section VI. C, contact law [25] that allows 
for some bounded tensile forces in the contacts, is such 



FIG. 24: Force ratio Q versus coordinate x of node 1 for the 
system of fig. 1221 without the ASD. 




FIG. 25: Same as figs. l9land [23l in the case of the system of 
fig. 1101 without the ASD. The straight lines are the tangents 
to the boundary curve at points A and B. 



that both equilibrium positions A and B will be simul- 
taneously possible for some values of Q, in the system of 
fig. [TQJ Q then varies with x exactly as shown on fig. [5H1 

with Qa = -p + — and Qb = 4= - — ■ 

V3 Fy ^3 Fy 

One may note, however, that the plasticity due to co- 
hesion of finite strength differs from the one due to geo- 
metric rearrangements in the two following respects. 
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FIG. 26: Same as figs. [2JJ and 1241 in the case of the system 
of fig. [10] outside the ASD. The force-displacement relation is 
now history-dependent, as shown by the arrows. 



• With contact law [25] plasticity does not disappear 
in the limit of small motions (when the ASD be- 
comes valid). 

• It is sensitive to the magnitude of external forces, 
not only on their direction. The figure analogous 
to 1261 in the (Q = x) plane, now depends on 
the value of F y . When F y is very much larger than 
/o, the cohesive strength of contacts might be ne- 
glected, and vanishes as a source of plastic dissipa- 
tion. 

It might be expected, on going, from the elementary 
example dealt with in this section, to larger and larger 
systems, that curves like fig. US] forces (like F) averag- 
ing to stresses and displacements (like x) to strains, will 
gradually look like fig. [20j In larger systems, the curve 
of fig. US] w iU look like a staircase. Presumably, as the 
system size increases, the number of the steps, and their 
amplitude, if expressed in terms of intensive quantities, 
will tend to zero. Then the smoothness of the curves 
sketched on fig. [20] might be recovered in the thermody- 
namic limit. Whether it actually will is of course not 
obvious a-priori, a careful statistical analysis f49j is re- 
quired. In the case of systems treated within the ASD, 
each step of the resulting staircase will be retraced back 
and forth, without any irreversibility. Such models can 
be expected to share the properties of the lattice system 
of ref. jl7| and paragraph VII. B. 8, in which the stair- 
case does indeed approach a smooth stress-strain curve 
in the thermodynamic limit. (But this curve is unique, 
one cannot obtain fig. [20] in such a case) . 



The difference between plasticity of cohesive and non- 
cohesive grains that was pointed out above is reminiscent 
of the difference in the behaviour, under growing hydro- 
static pressure, of sands and clays [jgj . As the magnitude 
of the load increases (but its direction is fixed), the level 
of plastic deformation in the cohesive material (clay) is 
much higher than in the non-cohesive one (sand). 

It is also interesting to note that some theories of fric- 
tion between solid surfaces [H[ are, just like the mech- 
anisms for internal friction that we invoke here, based 
on the history-dependent selection of one among several 
possible stable equilibrium configurations. 



C. Consequences of isostaticity. 

We focus here on systems of frictionless, cohesionless 
and rigid spheres (the contact law being [21]) in equi- 
librium under a given load, for which it was shown, in 
two steps (sections VI and VIII), that the force-carrying 
backbone is an isostatic structure. We discuss some spe- 
cific consequences of this property. In the simple exam- 
ple treated in subsection B just above, both equilibrium 
configurations A and B correspond to isostatic contact 
structures, and it is easy to predict for which value of 
the loading parameters the system will change from one 
to the other. Exploiting the isostaticity property, we will 
show here that such a prediction can, to some extent, be 
done in an arbitrary system. 

In this subsection, we only consider the backbone, ig- 
noring the rest of the system. We suppose that grains 
have been renumbered, so that index /i, with 1 < fi < Nf 
only label the degrees of freedom of objects that belong 
to the backbone. We shall also adopt the convention 
that the whole backbone does not move as a rigid body 
(thus excluding the fco corresponding degrees of freedom 
from the list). Likewise, 1 < I < N here only labels the 
force-carrying contacts (N — Nf). 

1. Response to perturbations, without rearrangement. 

Isostaticity of the whole structure means that matrix 
G, and its tranpose G T are square and have an inverse. 
Not only are equilibrium forces, given the load, uniquely 
determined, but it is also possible to predict how small 
external force increments (on the backbone) will be dis- 
tributed in the existing contacts. Changing the load from 
(F^ xt )i<^<N } to (F^ xt +SF^ t ) 1 <^< Nf will result, in con- 
tact I, in force increment Sfi, given by (summation over 
repeated indices implied) 

6ft = (G T )^6F^ = G-JSF^. (49) 

The backbone being rigid, this change in forces does 
not entail any displacement: = for each /i. This cor- 
rectly describes the mechanical response of the granular 
assemblage as long as all contacts forces remain positive. 
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This should be the case, in a finite system, for sufficiently 
small perturbations of the initial load. 

2. Dual response of velocities to bond length variations. 

Parallel to the one-to-one correspondence between con- 
tact forces and external loads expressed by eqn. is the 
inversible linear mapping between velocities and relative 
normal velocities in the contacts. There is no compati- 
bility condition in the absence of hyperstaticity, and one 
may impose arbitrary values to relative normal velocities 
(5Vi)i<i<n for the whole list of contacts. The result- 
ing velocities of the spheres are then (summation over I 
implied) : 

V» = G-l5V h (50) 

On comparing to I49[ it appears that the same matrix 
element G^ 1 is both equal to the force increment in con- 
tact I created when a unit external force is exerted on the 
coordinate \i on the one hand, and to the velocity coordi- 
nate /i when Sv is equal to one in contact I and to zero in 
all other contacts, on the other hand. Such a symm etry 
in response functions was remarked by Moukarzel [34| . 
who derived it by different means. 

3. Response to perturbations: structural rearrangements. 

The particular form of mechanical response expressed 
by eqn. 1491 in which no motion occurs and the load in- 
crement is supported by the initially existing contacts, 
ceases to be relevant as soon as negative contact forces 
appear. In the case of a two-parameter loading mode, 
such as the biaxial experiment at constant p, in which 
q is gradually increased from its initial value q = 0, one 
may write in each contact I 

fi = PiP + nq, 

where /3; and 7; are, due to isostaticity, geometrically 
defined coefficients. In general one finds that some of 
the 7; are negative. Let us denote as L~ the set of such 
contacts. The load will no longer be supported as soon 
as q reaches the value 

q m ax = mm p. (51) 

leL- 7; 

For larger g's, the theorem of virtual power shows that 
it is possible to decrease the potential energy upon open- 
ing the contact Iq for which the minimum in the right- 
hand-side of ED is reached, all other contacts remaining 
closed. The system will then rearrange, until a new set 
of contacts is created, such that the new load (p, q) is 
supported with positive contact forces. If one uses the 
ASD to describe this motion, then, within this approxi- 
mation, the new list of contacts, as shown in section VI, 



is entirely determined by the sole system geometry, as 
the solution to a simplex problem. Outside the ASD, the 
new equilibrium state, after the system rearranges, might 
depend on specific dynamical laws. In general, the range 
of validity of the ASD and the influence of the dynamics 
are to be tested, in experiments or, perhaps more easily, 
in numerical simulations. However, we have just shown, 
in fact, that the direction of velocities at the beginning of 
the rearrangement is determined by purely geometrical 
conditions, at least if ^ is unique: to find those direc- 
tions, just impose 5Vi = —1 (thus opening contact Iq), 
and SVi — for any I ^ Iq, from which all velocity com- 
ponents are deduced as w M = — G~ ; *, from equation [50] 

Simulations of disordered systems of discs [5(3] suggest 
that is generically unique, except in situations when 
the opening contacts involve a cluster of d+l-coordinated 
spheres in d dimensions. Examples of such clusters are 
sets of discs 8, 19 and 2, or 6 and 15, or 12 alone on 
figure [3] It is easily realized that once one contact force 
involving e.g., disc 8 is known, then all contact forces 
involving discs 8, 19, or 2 are also known, and propor- 
tional to the first one. Thus, they all vanish simultane- 
ously. This means that all matrix columns (G~l)i<fj,<N f 
are proportional to one another for all indices l that la- 
bel contacts of d-spheres belonging to the same d+l- 
coordinated cluster. Returning to the determination of 
the motion when the load ceases to be supported by the 
initial list of contacts, it follows that even though, in 
such a case, several contacts, involving the same cluster 
of d+l-coordinated spheres, may simultaneously open, 
the uniqueness of the initial velocities, up to a common 
amplitude factor, is preserved for all spheres that do not 
belong to the said cluster. 



4- Fragility. 

When a rearrangement occurs after a load increment, 
the mechanical response of the granular assembly, unlike 
the one expressed by equation 1491 involves both force 
changes and displacements. It depends on the possibil- 
ity of closing contacts that are not present in the initial 
equilibrium configuration. This geometric information is 
not contained in matrix G, which only depends on the 
network of initially existing contacts. One could thus 
study a second type of response to perturbations, that 
involves displacements. To see which of the two kinds of 
response is more relevant for the macroscopic mechanical 
behaviour, one has to impose perturbations that possess 
some macroscopic meaning, such as changes of q in a bi- 
axial experiment. Then, assuming, to fix notations, q is 
increased from zero, two cases need be considered. Either 
the thermodynamic limit ofq maxl as defined in 1 5 H is pos- 
itive, or it is equal to zero. In the first case, there exists a 
finite interval of stress for which no motion occurs in the 
continuum limit, and the mechanical response discussed 
in the preceding paragraphs in terms of the sole matrix 
G is macroscopically relevant. In the second case the 
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granular material might be appropriately termed frag- 
ile, since, in the thermodynamic limit, arbitrarily small 
macroscopic perturbations provoke rearrangements of the 
contact structure. Then, any macroscopic mechanical 
experiment involves displacements, the sole knowledge of 
one network of contacts that corresponds to a given value 
of the loading parameters is not sufficient. The response 
expressed by the sole matrix G is not the macroscopically 
relevant one. 

Our simulations of frictionless rigid discs [13, E, 
l5(ij show that such systems are indeed fragile in this 
sense. 
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5. An algorithm to compute a sequence of equilibrium 
configurations. 

This suggests the following procedure to determine the 
sequence of equilibrium states reached by an assembly 
of rigid, frictionless, cohesionless spheres under varying 
load (p,q), without resorting to any dynamical param- 
eter (without introducing any inertia, or mechanism of 
dissipation). 

• 1) Starting from an equilibrium configuration, in- 
crease loading parameter q until contact force fi 
vanishes. 

• 2) Move grains in the direction determined by the 
opening of contact Iq, the others remaining closed. 
Keep the same prescription for the grain trajecto- 
ries as for the initial velocities, taking into account 
the rotation of vectors n^ , until some new con- 
tact l\ is created, such that the new contact list, 
replacing Iq (now open) by l±, defines an isostatic 
structure. 

• 3) If, in the new contact structure, the contact 
forces that balance the load are all positive, a new 
equilibrium state, corresponding to the new load, 
has been reached: one may go back to step 1) and 
further increase q. Otherwise, some contact forces 
are negative. Pick up the one with the highest ten- 
sile force, call it Iq and go back to step 2), with the 
new contact list. 

This algorithm has been implemented by G. Combe 
and the present author We propose to name it the 
'geometric quasi-static method' (GQSM). It does involve 
arbitrary ingredients: there is no reason to forbid other 
openings of contacts once interstice hi has reached a 
finite positive value. Its great advantage is the possibility 
to compute trajectories from the sole knowledge of the 
system geometry. 

The system evolution, under a varying load, appears 
as a sequence of equilibrium states that are separated 
by 'jumps' or rearrangements, in which the list of active 
contacts is altered. In a phase of equilibrium, the forces 
are carried by a minimum list of contacts. In a phase 



of motion, normal relative velocities, among the whole 
bond list, are localized on one bond (several if a structure 
- a list of bonds- larger than the contact structure, is 
considered). Both maximum localization phenomena are 
related to geometric constraints. 

The predictions of the GQSM algorithm were com- 
pared with those of other methods that resort to dy- 
namical models (and, as argued in the introduction, also 
involve arbitrary, non-physical features). The results will 
be presented elsewhere. As mentioned above, mechanical 
properties, at the level of individual trajectories in con- 
figuration space, cannot be expected, outside the ASD, to 
be uniquely determined. However, in view of the impor- 
tant role of the geometry, which determines exactly the 
value of the loading parameters for which system should 
rearrange and the direction of the initial velocity vector, 
it can be hoped that the statistical properties of such tra- 
jectories that are relevant for the macroscopic laws will 
present little dependence on dynamical features of the 
system (such as masses or dissipative shock laws) . 



6. Rearrangements within the ASD. 

Within the approximation, as the equilibrium state 
corresponding to a given load is unique, there is no need 
to resort to an incremental approach. If one however does 
so, then the whole rearrangement event is geometrically 
determined. It can be computed with the GQSM as pre- 
sented above. Then, it will be observed, on performing 
step 3) of the algorithm, that the new contact structure, 
as soon as a new contact is created, supports the load 
with only positive contact forces. Thus, unlike in the 
general case [50( , no cascade of successive rearrangements 
occurs in step 3). Rearrangements are simpler events in 
which one element of the contact structure is replaced by 
another. 

Let us prove this statement. 

Let So denote the old list of contacts, and Si the new 
one. Both structures are isostatic, and for any given load 
one can find unique values of both sets of bond forces 
(fi)ieS an d (fi)ieSi that ensure equilibrium. In the fol- 
lowing members of these two sets, in order to distinguish 

them, are written down with a superscript: /, and //^ 
respectively denote the force carried by bond I, as com- 
puted with structure Sq and with S\. 

Recalling also the notations of the preceding para- 
graph, Si is equal to So, deprived of contact Iq, to 
which contact li is added. When the value q max of 
the loading parameter is reached, f® has decreased to 
zero. This means that, exceptionally, the smaller struc- 
ture So \ {Iq} = Si \ {h} can support the load, and one 



has while 



// 0) for each I € S \ {k}- As 



we assume, for simplicity, that contact forces reach zero 
separately, there exists a finite range of positive incre- 
ments Sq such that one has / ; < 0, while /j > for 
I G Sq \ {Iq}, for q = q max + Sq. Likewise, reducing the 
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> 



Sq interval if needed, we require the condition /; 
for leSi\ {k}. 

We now pick up one such value of q, and evaluate the 
variation SW of the potential energy (that corresponds 
to this value of q) in the rearrangement. 

On the one hand, one may obtain SW on applying the 
theorem of virtual work to structure So- As contact Iq has 
opened, the corresponding relative normal displacement 
is negative: Sui a < 0, while Sui — for each I =/= l a . 

Therefore, because f, < 0, one has 

SW = -f^Su l0 < 0. 

On the other hand, one may obtain SW on applying the 
theorem of virtual work to structure Si. As contact li has 
closed, the corresponding relative normal displacement 
is positive: 8ui t < 0, while Sui — for each / =/= li. 

Therefore, because SW = —fj/Suu < 0, one has 



/?> > o. 

Thus, the new contact structure supports the load with 
positive contact forces as soon as q > q ma x, and a new 
stable equilibrium state has been reached. 

In the general case, we stressed the difference between 
the mechanical response of the granular system without 
rearrangement, which can be deduced from the geometry 
of the contact structure, via matrix G, and the mechan- 
ical response involving some rearrangement, the deter- 
mination of which requires some additional prescription 
(such as that of the GQSM) to move the particles. 

This difference is much less important within the ASD: 
as the matrices G pertaining to either structure do not 
change in the motion, all displacement coordinates will 
simply be found as follows: 



= h 



(52) 



where h\ x denotes the initial opening of contact l\ and 
the matrix G is that of structure Si . 

Equations [50] and [55] only differ by a scale factor, in- 
terstice hi x . There is nothing especially singular in the 
distribution of open interstices in dense granular systems 
at equilibrium. So, it can be expected that macroscopic 
averages corresponding to both response functions, \EU\ 
and 1521 are proportional to one another. Moreover, the 
response without rearrangement, expressed bv !50l is the 
same with and without the ASD. 

In the following subsection, we derive explicitly the 
form of the macroscopic response function to small in- 
crements in applied external forces, in the case of the 
triangular lattice model. These are large scale averages 
of (combinations of) microscopic responses expressed by 
eqn.[H 

We shall therefore speculate that the results to be de- 
rived below, for the form of such macroscopic Green's 
functions, are also valid for the average of response func- 
tions without rearrangements in general. 



D. Macroscopic response of the triangular lattice 
model. 

In the model system studied in ref. [l7|, the results of 
which are recalled in paragraph VII. B. 8, it is possible to 
find the form of macroscopic equations to be solved when 
a small density of external forces Sf ext is superimposed 
over an initial equilibrium state. 

To do so, one just needs to translate the properties 
stated in paragraph VII. B. 8 in incremental form. 

First, let us impose, without loss of generality, a few 
conditions on function / defined in !45l It is convenient to 
choose a symmetric function of e a p and ep a , the deriva- 
tion in [35] being taken regarding both strain components 
as independent variables. 

Then, defining, in e space, a norm ||e|| by 

H||| 2 =|:|=eii+2ei2 + 4, 
one may enforce (replacing / by //||V/||) the condition: 



|V/|| = 1, 



(53) 



everywhere on E. 

One starts from an equilibrium state in which the stress 
field, a, is assumed to stay strictly inside the supported 
range, defined by inequalities 1421 everywhere in the sys- 
tem. This initial state is also characterized by a displace- 
ment field Uo and a strain tensor field e (everywhere on 
E, and abiding by [46]), the origin being defined by the 
reference state (the undisturbed regular lattice of spac- 
ing a). One then looks for the stress increment field §a_, 
displacement increment field u and strain increment field 
St that result from the application of Si ext . The problem 
is dealt with to first order in any of these quantities, that 
are linear in 8f ext , assumed small. 

Let us define 
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dta/sde^s ' 



a fourth-order tensor that depends on e. One has, upon 
differentiating the macroscopic law 



A 



df 



de, 



a/3 



the decomposition of stress increments as 

6<r aP = Su^l +8a ( V, 
with (summation over repeated indices) 

8&ap = ^Aa^sSe-yg, 

and 5<t%£ = tWs- 

Condition [53] yields, by derivation, 

Of 



de a f 



0, 
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whence the orthogonality between a and Sa^. Since e 
must remain on E, 5e is also orthogonal to a. 

In view of the symmetry of the stress tensor and of the 
conditions imposed on function /, tensor A satisfies the 
following symmetries: 

Because it is a second-order derivative, one also has: 

A a f3yS = AySafl- 

Tensor A is thus endowed with the same symmetry prop- 
erties as a tensor of elastic constants (or of viscosity co- 
efficients). 

We have seen that it might be viewed as a linear oper- 
ator within the space of symmetric second-order tensors 
that are orthogonal to a, or, in other words, within the 
tangent plane to surface X in strain space. Because of 
the strict convexity of T>, this operator is positive definite 
(this is easily realized, as the curvature of £ is turned in- 
wards). 

Transforming the equilibrium equation into one for the 
unknowns u and <5A, using 143] one obtains (d a denoting 
a derivative with respect to coordinate a) 

dp [XA a ^ s d s u-y] - dp (^trj\ + &fl xt = 0, (54) 

while the displacement field should satisfy 

(Japdpua = 0. (55) 

Equations [54ll55l supplemented by suitable boundary 
conditions, define, because of the positive-definiteness of 
operator A, an elliptic boundary value problem. The so- 
lution is unique provided 2 conditions (in 2D) involving 
u and/or its normal derivatives are specified everywhere 
on the system boundary. 

We now turn to the situation when the initial stress 
field is a uniform hydrostatic pressure: 

with a position-independent pressure Pq. In view of con- 
dition [53] on /, it should be noted that A coincides with 
PoV% in this case. The corresponding tangent space to 
£ is the space of traceless tensors. 

In general, tensor A reflects the common symmetries of 
the material (the triangular lattice) and the stress tensor. 
In this particular case, it will possess all the symmetries 
of the regular triangular lattice. The tensor of elastic 
constants, in that case [511 ] . has the same symmetries as 
in an isotropic medium. Because it operates within the 
space of traceless tensors, tensor A reduces to a scalar 
K: one has, for any traceless strain increment, 

Aaf3~/80~£~/S — Kde a p. 

[SUhas become 

K PoV2V 2 u - V(5P) + 5f ext = 0, 



while [55] now states that the displacement field should be 
divergenceless: 

V ■ u = 0. 

One recognizes the Stokes problem for viscous incom- 
pressible flow, in which the displacement replaces the 
velocity field, the product KPq^/2 plays the role of the 
shear viscosity, and SP is a pressure field to be deter- 
mined on solving the full boundary value problem. 

Green's functions for the Stokes problem can be found, 
e.g., in [52l |. In an infinite 2D medium, the velocity 
field varies logarithmically with the distance to the point 
where a concentrated force is applied. 



E. Discussion. 

From the results just above, it can be concluded that 
the form of the macroscopic equations ruling the dis- 
placement field created by a small perturbation to a pre- 
stressed granular sample in equilibrium should be elliptic, 
provided the microscopic rearrangements are dealt with 
within the ASD. 

From the discussion at the end of paragraph IX. C. 6, we 
expect that operator G -1 , in the general case, also aver- 
ages macroscopically as the Green function of an elliptic 
second-order partial-differential operator. One may ob- 
tain a suitable macroscopic average on taking, e.g., the 
mean of all matrix elements G^ 1 for which the vector 
pointing from bond I to the center of the grain which co- 
ordinate [i belongs to is in some prescribed small neigh- 
bourhood of a given vector. 

G -1 rules the response without rearrangement. The 
general -and, in view of the fragility property, most 
relevant- case of mechanical response involving rear- 
rangements outside the ASD appears to involve more ge- 
ometric information than the one contained in matrix G: 
it could be observed [50( that step 3) of the GQSM al- 
gorithm introduced in paragraph IX. C. 5 could involve 
a long sequence of elementary rearrangements replac- 
ing one contact by another. Unlike the distribution of 
open gaps between adjacent particles, that of the magni- 
tude of such complex rearrangements can be quite wide 
and might significantly affect the macroscopic response 
in terms of dispacements. This will be studied in a forth- 
coming publication. In the case of a disordered granular 
assembly, no small parameter, like the level of polydis- 
persity of discs in the triangular lattice model, is avail- 
able to control the validity of the ASD. As found in sec- 
tion VIII, stable equilibrium states of frictionless discs 
or spheres are especially scarce in configuration space, as 
full rigidity is required. Outside the ASD, impenetrabil- 
ity constraints do not limit a convex accessible domain of 
configuration space. Whereas the route from one equilib- 
rium state to another, within the ASD, can be straight, 
it might have to follow a long and tortuous path out- 
side the approximation. (The ASD amounts to simplify 
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this complex geometry, straightening up local curvatures, 
etc. . . ) 

Interestingly, Tkachenko and Witten [35(, following a 
suggestion by Alexander [4(|, speculated that, as a conse- 
quence of the isostaticity property, the mechanics of fric- 
tionless sphere packings should be described, at the con- 
tinuum level, by laws of the type proposed in refs. [H>[l3|: 
the response to perturbating force fields satisfies hyper- 
bolic partial differential equations. From considerations 
on the floppy modes that appear within a subsystem that 
is isolated from the rest of the sample, they derive a simi- 
lar directional structure for matrix G _1 as for the macro- 
scopic response in such theories: in a sample limited by a 
free surface in the upwards direction, force perturbations 
are not felt above the point where they are introduced. 

Although we do not venture here to speculate on the 
form of macroscopic equations that rule the mechanical 
response with rearrangements in a general, disordered 
system for which the ASD might not be valid, our con- 
clusions above do go far enough as to clearly contradict 
the ones of [35] , since those are concerned with the same 
object (operator G _1 ). 

An explanation for this discrepancy could be that 
Tkachenko and Witten mainly based their conclusions 
on the observation of packings (numerically) obtained by 
sequential deposition algorithms under gravity. 

When the stress tensor approaches the boundary of 
the region of supported loads (i.e., when one of the con- 
ditions in 221 is almost an equality) one can observe [l7l ]. 
for the triangular lattice model, that the list of force- 
carrying contacts approaches a limit that comprises all 
the bonds parallel to two of the three lattice directions, 
and none of the bonds parallel to the third. The topology 
of the backbone thus approaches that of a square lattice. 
In this particular case [35], it is easy to check that a 
description in terms of force propagation, involving hy- 
perbolic equations, applies. The marginally supported 
stress states of this model are the analog of the Coulomb 
condition for an isotropic medium. When the Coulomb 
criterion is everywhere satisfied as an equality, the ma- 
terial is everywhere on the verge of plastic failure, and 
it has long been known (and exploited for the evaluation 
of critical loads [53]) that the macroscopic equations are 
of the hyperbolic type. This situation has been termed 
'incipient failure everywhere' (IFE) in [H, [Hj]. 

One may conjecture that deposition algorithms [5^. [55| 
will systematically produce internal states close to IFE. 
Specifically, we expect sequential deposition under grav- 
ity to result in the 'active' Rankine state, in which the 
pressure on the lateral walls is barely sufficient to contain 
macroscopic plastic flow of a horizontal granular layer 
sumitted to its own weight. In the case of discs with a 
small or moderate polydispersity in 2D, the deposition 
algorithms do in fact produce networks of force-carrying 
contacts that are very close to the limiting states of the 
triangular lattice model (a deformed square lattice). 

Therefore, we suspect that Tkachenko and Witten's 
arguments only apply to those particular cases of limit 



states or IFE. 

There are, apart from the arguments put forward 
in [35j, other aspects on which the general properties 
we have been discussing as well as the numerical re- 
sults obtained on the triangular lattice model appear at 
odds with the assumption of a direct relationship between 
stress components, and related theories. Leaving a more 
complete discussion to subsequent work, let us merely 
point out that the nature of the boundary conditions has 
dramatic effects if the macroscopic equations are hyper- 
bolic. In fact, if a rigid boundary transmitting a stress 
is replaced by a distribution of external forces imposed 
independently on the grains that are close to the edge, 
such theories predict this change to significantly affect 
the whole system (which has lost its rigidity). In our ex- 
perience 0, [IB], some rearrangement does occur, but its 
effects arc confined to a boundary layer of finite depth. 

We also note that our results disagree with some of 
Moukarzel's [H, [HI, predicting perturbations due to a 
localized force to increase exponentially with distance. 
Although his results are very accurate and were obtained 
on very large systems, the propagative nature of forces, 
which can be calculated from 'top' to 'bottom' in a sin- 
gle sweep, is an explicit ingredient of his model, that 
was adapted from the one of [3^]. Our results on the 
triangular lattice model disagree with his because this 
very large effect of force perturbations (or, equivalently 
- see 032 an d I5"D1 - of bond length variations) would cause 
the level of distortion of the regular lattice, due to the 
polydispersity of discs, to increase very fast with the sys- 
tem size. Rather, we observed it to approach a finite 
thermodynamic limit. Once again, we suspect that the 
very peculiar properties obtained in these studies stem 
from the consideration of a special case in which forces 
happen to possess a propagative nature. 

Finally, the (provisional) conclusion we propose here is, 
as already mentioned in section VIIE, that the rigidity 
of the grains and the isostaticity property do not neces- 
sarily entail very special, critical or singular macroscopic 
mechanical properties. Moreover, we expect - as systems 
dealt with within the ASD exhibit the same kind of elas- 
ticity as networks of rigid cables - that if unusual, exotic 
properties exist, then they are related to the displace- 
ments (the rearrangements) rather than the network of 
forces (or the operator G attached to it). 



X. CONCLUSION AND PERSPECTIVES. 

Let us first briefly summarize the main results pre- 
sented in this paper. 

Specializing to frictionless grains, and assuming that 
granular packings, under slowly varying sollicitations, 
tend to stable equilibrium states, we have shown that 
geometry determines, to a large extent, the mechanical 
behaviour of such materials. 

Spatial arrangements of granular packings in equilib- 
rium under a given load are quite specific points in config- 



38 



uration space. Rigid grains that only exert normal con- 
tact forces on one another, once submitted to a supported 
load, will generically pack in such a way that the prob- 
lem is isostatic, i.e., there is no indeterminacy of forces. 
The value of all contact forces is determined by equilib- 
rium equations and the geometry of the contact struc- 
ture. This yields a rigourous upper bound on the con- 
tact coordination number of any packing of rigid grains. 
These properties hold for compressive or tensile contact 
forces. Contact structures, in equilibrium, are not al- 
ways rigid, especially (but not exclusively) in the case 
when contacts can sustain tensions. Even if loose parti- 
cles, that carry no force, are discarded from the count, 
the upper bound on the coordination number might not 
be reached. 

If the packing is such that the approximation of small 
displacements might be well justified, in particular in 
the case of regular arrangements on lattices, stronger 
properties were established, provided the problem can be 
coped with in the framework of convex optimization the- 
ory (which requires the definition of a potential energy, 
thus excluding finite strength cohesion). Then 

• Not only the forces once the contact structure is 
known, but the force-carrying structure itself is en- 
tirely determined by the system geometry. 

• Grain positions are also determined, apart from 
possible 'floppy mode' motions, of bounded ampli- 
tude, that do not affect the value of the potential 
energy. 

• Displacements from the reference configuration on 
the one hand, and contact forces on the other hand 
are the solutions to two optimization problems in 
duality. 

• For rigid grains, force-carrying structures are the 
exact analog of cost-minimizing directed paths in 
scalar transport problems. 

Such situations are thus very attractive from a theorist's 
point of view: the reduction of the mechanical problem to 
one of random geometry is complete, and analogies with 
other models of theoretical statistical physics (directed 
percolation, directed polymer in a random environment) 
can be drawn and exploited. However some important 
features of granular mechanics are absent: such systems 
are devoid of plasticity and hysteresis. 

Pursuing the stability analysis beyond the ASD in the 
case of discs or spheres, we have shown that the force- 
carrying structure must be rigid if contacts do not with- 
stand tension, because any floppy mode would imply in- 
stability. This entails that the force-carrying backbone in 
systems of rigid spheres is, generically, an isostatic struc- 
ture, its coordination number is equal to 2d in dimension 
d. 

Analogous systems of cables (that resist tension, but 
no compression), on the other hand, will generally keep 



some amount of Soppiness, since mechanisms in the equi- 
librium state are all stable. 

Assemblies of frictionless grains will, in general, exhibit 
internal friction, due to the multiplicity of stable equi- 
librium states corresponding to the same external load. 
This non-uniqueness might stem from the finite extent of 
rearrangements or from bounded cohesion forces. 

If submitted to slowly varying loads, packings of rigid 
grains will evolve via a succession of jumps or crises sepa- 
rated by phases of rest. The isostaticity property implies, 
for a system of rigid frictionless spheres, that the concen- 
tration of forces is maximal during a phase of rest (forces 
cannot be carried by a strictly smaller set of contacts), 
and that the concentration of deformation is maximal at 
the beginning of a jump (there cannot exist a strictly 
smaller list of interstices in which relative normal veloc- 
ities are not equal to zero). 

Although the motion in a rearranging event depends on 
the actual granular dynamics, the forces during a phase 
of rest, and the direction of velocities at the beginning of 
motion, are geometrically determined. 

Two kinds of response functions to force increments 
can be studied, depending on whether the perturba- 
tion provokes a change in the contact list. Some recent 
studies of response functions, without rearrangement of 
the grains, were discussed and we argued that some of 
their conclusions might be specific to sequential deposi- 
tion models, in which forces can be propagated along a 
preferred direction. The fragility of frictionless granular 
assemblies in the thermodynamic limit implies however 
that macroscopically meaningful perturbations always in- 
volve some amount of rearrangement. 

The results of the present article suggest both general 
perspectives and specific problems, to be dealt with in 
future work. 

An important feature of granular materials is the spar- 
sity, in configuration space, of equilibrium configurations. 
Those, especially for rigid grains, have very specific char- 
acteristics. Moreover, they are generally suitable for one 
particular load. In such circumstances, it might not be 
adequate to choose first one specific geometric arrange- 
ment and contact structure, built, e.g., by some conve- 
nient algorithm that respects impenetrability conditions, 
and then to apply external forces and see how they could 
be balanced by contact forces. The list of active contacts 
is itself chosen according to the external load. Many re- 
cent studies were devoted to the way forces distribute 
among a fixed list of contacts, and to the ensuing statis- 
tics of contact force values. Although models along these 
lines might capture some of the physics, they ignore dis- 
placements. Displacements, as our results have amply 
shown here, are always part of the problem. The very 
definition of a force requires the consideration of some 
amount of displacement. A normal reaction force in the 
frictionless contact between two rigid objects is a geomet- 
rically defined quantity, a Lagrange parameter associated 
with an impenetrability constraint in configuration space. 
Large assemblies of frictionless rigid grains are fragile: 
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tiny load increments will be associated with rearrange- 
ments of the contact structure. If one wishes to under- 
stand the macroscopic mechanical behaviour of granular 
systems and its relationship to grain-scale phenomena, 
the question of the magnitude of such rearrangements, 
in which the system moves from an equilibrium state to 
another, is crucial. 

Other, more specific questions, that are related to 
statistics and the continuum limit, naturally follow from 
the mechanical properties we have been presenting. 
When is the ASD is a good approximation, apart from 
lattice models ? Are the same states periodically revis- 
ited in cyclic sollicitations ? What will be the density 
and the effect of floppy modes in systems of non-spherical 
frictionless particles ? Will the staircase-like stress-strain 
curve approach a smooth limit when the system size in- 
creases ? To what extent are rearrangements sensitive to 
the actual dynamical rule ? Such problems would benefit 
from careful numerical simulations, and we shall address 
some of these questions in forthcoming publications. 

The treatment of granular systems with friction could 
be tackled with a similar approach to the one developped 
here: one could investigate the range of stability of a 



given contact structure, as the load gradually varies, by 
purely static means. In the presence of friction, granular 
packings are also observed, in experiments and dynamic 
numerical simulations, to evolve by a succession of crises 
localized in time. We expect the geometry of the assem- 
blage to dictate, to a large extent, the way such sudden 
motions are initiated. 

It can be concluded that much of the promising 
prospects, as well as much of the difficulties ahead, in 
the study of mechanical properties of granular materi- 
als close to equilibrium, are in the understanding of the 
disordered, yet quite peculiar, geometry of large systems 
that adapt their contact network to sustain the load. 
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